微分和差分方程Word文档下载推荐.docx
《微分和差分方程Word文档下载推荐.docx》由会员分享,可在线阅读,更多相关《微分和差分方程Word文档下载推荐.docx(34页珍藏版)》请在冰豆网上搜索。
[x,fval]=fsolve(fun,x0)
其中,fun为函数名,x0为初始值,options为参数设置(迭代算法、步长、精度等),x为最优解,fval为最优解对应的最优值。
示例如下:
x0=[-5;
-5];
%随机初始值
options=optimset('
Display'
'
iter'
);
%显示每次迭代的过程
[x,fval]=fsolve(@myfun,x0,options)%调用优化函数
其中,myfun如下:
functionF=myfun(x)
F=[2*x
(1)-x
(2)-exp(-x
(1));
-x
(1)+2*x
(2)-exp(-x
(2))];
迭代过程:
最后的解为:
x=
0.5671
fval=
1.0e-006*
-0.4059
该函数求解性能很强,大家可以随意设置初始值x0看看结果如何,理论上结果是不变的,也就是该函数使用的迭代算法对初始值不敏感。
(4)符号方程——solve
solve(eq)
solve(eq,var)
solve(eq1,eq2,...,eqn)
g=solve(eq1,eq2,...,eqn,var1,var2,...,varn)
其中,eq为待求解的符号方程,默认的求解符号是x;
通过var设置带求解的符号变量;
求解符号方程组的话,需要依次添加每个方程,以及带求解的所有符号变量。
g是存放结果的结构体变量,通过g.x调用显示结果。
%%求解a*x^2+b*x+c=0关于x的符号解
solve('
a*x^2+b*x+c'
)
%%求解a*x^2+b*x+c=0关于b的符号解
b'
%%求解方程组x+y=1,x-11*y=5的所有符号解
S=solve('
x+y=1'
x-11*y=5'
S.x,S.y
%%求解方程组a*u^2+v^2,u-v=1,a^2-5*a+6=0的所有符号解
A=solve('
a*u^2+v^2'
'
u-v=1'
a^2-5*a+6'
A.a,A.u,A.v
结果就不写出来了,大家测试下就好了。
注意:
像这样a*x^2+b*x+c的写法代表方程a*x^2+b*x+c=0,所以u-v=1也可以写成u-v-1哦,不妨测试下。
(5)牛顿迭代法解非线性方程
牛顿迭代算法原理比较简单,大家可以参考Lagrange定理推导。
求解示例如下:
%%使用默认精度1e-6
x=newton('
myfun'
[1;
2],10)
%%设置精度
2],50,0.0000001)
其中newton函数如下:
functionx=newton(funname,x0,Maxgen,tol)
%作用:
通过牛顿迭代算法求解非线性方程组
%调用方式:
x=newton('
f_name'
x0)
%x=newton('
x0,tol)
%
%x:
最优解
%funname:
定义方程组的函数名
%x0:
初始值
%Maxgen:
最大迭代次数
%tol:
精度(默认:
1e-6)
h=0.0001;
M=Maxgen;
%最大迭代次数
ifnargin<
4
tol=1e-6;
end
x=x0;
xb=x-999;
n=0;
whileabs(x-xb)>
tol
xb=x;
ifn>
M
break;
end
y=feval(funname,x);
fprintf('
n=%3.0f,x=%12.5e,y=%12.5e,\n'
n,x,y)
y_driv=(feval(funname,x+h)-y)/h;
x=xb-y./y_driv;
n=n+1;
fprintf('
n=%3.0f,x=%12.5e,y=%12.5e,'
ifn>
\n'
disp('
Warning:
iterationsexceedsthelimitation,probabledevergent'
通过调节最大迭代次数很容易发现,一般在10次左右就可以达到最优解,或者近似最优解;
而通过加大迭代次数或者精度也很难获取更有解。
因而牛顿迭代可以通过很少的迭代次数获取较优解。
(6)遗传算法
这里讨论的遗传算法是基于matlab自带工具箱的,后面讲详细介绍GA算法的实现过程。
工具箱调用示例:
x=ga(fitnessfcn,nvars)
x=ga(fitnessfcn,nvars,A,b)
x=ga(fitnessfcn,nvars,A,b,Aeq,beq)
x=ga(fitnessfcn,nvars,A,b,Aeq,beq,LB,UB)
x=ga(fitnessfcn,nvars,A,b,Aeq,beq,LB,UB,nonlcon)
x=ga(fitnessfcn,nvars,A,b,Aeq,beq,LB,UB,nonlcon,options)
x=ga(problem)
[x,fval]=ga(...)
[x,fval,exitflag]=ga(...)
其中,fitnessfcn是自适应函数或者目标函数;
nvars是带求解变量的个数;
A和b分别是线性不等约束条件
的矩阵系数A和向量b;
Aeq和beq分别是线性等式约束条件
的矩阵系数Aeq和向量beq;
LB和UB分别为nvars个带求解变量的上下限向量。
这里面的向量都是列向量哦!
其他参数可以参考matlab的help。
应用举例:
目标函数:
约束条件:
求解程序:
%%线性非等式约束
A=[11;
-12;
21];
b=[2;
2;
3];
%%变量下限
lb=zeros(2,1);
%%调用工具箱
[x,fval,exitflag]=ga(@gatestfun,2,A,b,[],[],lb)
其中,gatestfun如下:
functiony=gatestfun(x)
p1=0.5;
p2=6.0;
y=p1*x
(1)^2+x
(2)^2-x
(1)*x
(2)-2*x
(1)-p2*x
(2);
结果如下:
0.66701.3340
-8.2258
exitflag=
1
exitflag是算法终止标志,1代表正常终止,其它参考help。
这里的函数参数传递是通过@gatestfun匿名方式的,而上面我们编写的newton是通过'
字符串方式传递的。
ga工具箱求解的是函数的最小值。
(7)黄金分割法和插值法——单变量
x=fminbnd(fun,x1,x2)
x=fminbnd(fun,x1,x2,options)
[x,fval]=fminbnd(...)
[x,fval,exitflag]=fminbnd(...)
[x,fval,exitflag,output]=fminbnd(...)
其中,fun是待求解方程的名称,x1和x2分别为最优解的上下限,options和前面一样参数配置,其他参数类似。
应用举例:
求解方程
在区间(0,2)上面的最小值,程序如下:
%%带求解函数的匿名表达
f=@(x)x.^3-2*x-5;
%%调用函数求解
x=fminbnd(f,0,2)
(8)梯度最优化算法——单、多变量
同样是求函数的最小值,调用示例:
x=fminunc(fun,x0)
x=fminunc(fun,x0,options)
[x,fval]=fminunc(...)
其中,fun为带求解函数名,x0为初始值,options为参数设置选项。
%%简单函数测试,不带设置
x0=[1,1];
[x1,fval]=fminunc(@myfun81,x0)
%%带设置的函数测试
options=optimset('
GradObj'
on'
[x2,fval]=fminunc(@myfun82,x0,options)
%%匿名函数测试
f=@(x)sin(x)+3;
x3=fminunc(f,4)
运行过程中可能出现warning信息,可以在程序中加入warningoff关闭警告信息。
(9)Nelder算法——多变量
x=fminsearch(fun,x0)
x=fminsearch(fun,x0,options)
[x,fval]=fminsearch(...)
[x,fval,exitflag]=fminsearch(...)
[x,fval,exitflag,output]=fminsearch(...)
其中,参数大家可以参考help。
%%匿名函数定义
banana=@(x)100*(x
(2)-x
(1)^2)^2+(1-x
(1))^2;
x0=[-1.2,1];
[x,fval]=fminsearch(banana,x0)
这几个函数调用方式几乎一样。
前面说到,初始值的选取对结果影响很大,即初始值选取是牛顿迭代算法的关键因素。
但是像遗传算法等智能算法对初始值参数并不敏感,甚至无所惧。
那大家应该有疑问了,既然如此,何必不适用遗传算法呢?
其实,如果深究这些算法的话,会发现:
对于大部分的方程求解,牛顿迭代最多只需要10次迭代就可以达到最优解,而其它算法既耗时也耗存储。
在实际应用中,我们既要考虑到精度又要兼顾计算复杂度,所以多数人宁愿在牛顿迭代初始参数选取上下功夫,也很少去研究如何改进其它算法性能,对于方程求解而言。
除此之外,还可以通过最速降、梯度下降、贪婪算法、模式搜索、爬山算法、遗传算法、神经网络等优化方法解决方程求根问题。
这些将在后续逐个讨论。
1.2常微分方程求解
参考附件:
《常微分方程的解法.pdf》、《常微分方程(ODEs)的MATLAB数值解法.pdf》,附件为常微分方程的求解方法和程序。
下面介绍两种常用的龙格库塔(Runge-Kutta)方法,实质上是求积分。
(一)定步长四阶龙格库塔法
1、算法原理
一阶常微分方程初值问题
(2.1)
的数值解法是近似计算中很重要的部分。
常微分方程初值问题的数值解法是求方程(6.1)的解在点列
上的近似值
,这里
是
到
的步长,一般略去下标记为
。
常微分方程初值问题的数值解法一般分为两大类:
(1)单步法:
这类方法在计算
时,只用到
、
和
,即前一步的值。
因此,在有了初值以后就可以逐步往下计算。
典型方法如龙格–库塔
方法。
(2)多步法:
时,除用到
以外,还要用
,即前面
步的值。
典型方法如Adams方法。
经典的
方法是一个四阶的方法,它的计算公式是:
(2.2)
方法的优点是:
单步法、精度高,计算过程便于改变步长,缺点是计算量较大,每前进一步需要计算四次函数值
2、程序实现
functionmain
%测试主函数
%判断测试方程存在与否
ifexist('
myfun.m'
)==0
没有为方程创建名为myfun.m的函数文件,请参照下例建立它'
functionz=myfun(x,y)'
z=y-2*x/y;
'
并将该文件保存在work文件夹下'
X1=input('
请输入求解区间的左端点X1='
Y1=input('
请输入微分方程的初始条件Y1=(X=X1时Y的值)'
Xn=input('
请输入求解区间的右端点Xn='
h=input('
请输入求解步长h='
X=X1;
Y=Y1;
%运算初始点
n=0;
%节点序号变量置零
whileX<
=Xn-h
K1=myfun(X,Y);
K2=myfun(X+h/2,Y+K1*h/2);
K3=myfun(X+h/2,Y+K2*h/2);
K4=myfun(X+h,Y+K3*h);
X=X+h;
Y=Y+h*(K1+2*K2+2*K3+K4)/6;
%四阶标准的龙格-库塔公式
n=n+1;
%节点序号加1
第%d个点的计算结果为X=%10.8f,Y=%10.8f\n'
n,X,Y);
plot(X,Y,'
o'
holdon
3、应用举例
测试的微分方程程序如下:
functionz=myfun(x,y)
%%微分方程为:
dy/dx=y-2*x/y
z=y-2*x/y;
4、测试结果
请输入求解区间的左端点X1=1
请输入微分方程的初始条件Y1=(X=X1时Y的值)3
请输入求解区间的右端点Xn=5
请输入求解步长h=.005
第1个点的计算结果为X=1.00500000,Y=3.01169404
第2个点的计算结果为X=1.01000000,Y=3.02344308
第3个点的计算结果为X=1.01500000,Y=3.03524747
......
(二)自适应变步长龙格库塔法
单从每一步看,步长越小,截断误差就越小,但随着步长的缩小,在一定求解范围内所要完成的步数就增加了;
步数的增加不但引起计算量的增大,而且可能导致舍入误差的严重积累。
因此同积分的数值计算一样,微分方程的数值解法也有个选择步长的问题。
在选择步长时,需要考虑两个问题:
1>
怎样衡量和检验计算结果的精度?
2>
如何依据所获得的精度处理步长?
考察经典的四阶龙格-库塔公式(2.2),从节点
出发,先以
为步长求出一个近似值
,由于公式的局部截断误差为
,故有:
(2.3)
然后将步长折半,即取
为步长从
跨两步到
,再求得一个近似值
,每跨一步的截断误差是
,因此有:
(2.4)
比较(2.3)式和(2.4)式我们们可以看到,步长折半后,误差大约减少到
即有:
由此得到下列的估算公式:
这样,可以通过检查步长,折半前后两次计算结果的偏差:
来判定所选的步长是否合适。
具体的说,将区分以下两种情况处理:
1、对于给定的精度
,如果
,反复将步长折半进行计算,直至
为止;
这时取最终得到的
作为结果。
2、如果
,反复将步长加倍,直到
为止,这时再将步长折半一次,就得到所要的结果。
这种通过加倍或折半处理步长的方法称为变步长方法。
self.m'
把work文件夹里没有self.m文件'
eps=10^(-8);
x1=input('
请输入初始点x1='
y1=input('
请输入初始条件y1='
xn=input('
请输入终止条件xn='
h1=input('
请输入初始步长h1='
h=h1;
h=%10.8f,x=%10.8f,y=%10.8f\n'
h,x1,y1);
%输出初始条件
ifh>
abs(x1-xn)
初始步长取得过大,超过了求解区间的长度'
x1=xn+1;
%终止运算
whilex1<
=xn
[u2,v2,h,err]=self(x1,y1,h);
H=h;
half_err=err;
double_err=err;
%当误差过大时,反复缩小步长
whilehalf_err>
eps
[u2,v2,h,err]=self(x1,y1,H);
%当误差过小时,不断将步长增大
whiledouble_err<
H=2*H;
%误差合适,最后进行调整运算
ifdouble_err>
=eps
H=H/2;
%输出此点结果,为下一节点的计算提供初始值
H,u2,v2);
x1=u2;
y1=v2;
h=h1;
plot(x1,y1,'
其中,self.m如下:
function[u2,v2,h,err]=self(x1,y1,h)
%
%该函数用来调整自适应
u1=x1;
%u1为x1的备份,供步长为h/2时计算下一个节点时使用
v1=y1;
%v1为y1的备份,供步长为h/2时计算下一节点数值解时使用
%用四阶经典公式计算步长为h时第1个节点处的数值解
k1=myfun(x1,y1);
k2=myfun(x1+h/2,y1+h*k1/2);
k3=myfun(x1+h/2,y1+h*k2/2);
k4=myfun(x1+h,y1+h*k3);
y2=y1+h*(k1+2*k2+2*k3+k4)/6;
%四阶经典公式计算步长为h/2时的第一个节点处的数值解
h=h/2;
fori=1:
2
k1=myfun(u1,v1);
k2=myfun(u1+h/2,v1+h*k1/2);
k3=myfun(u1+h/2,v1+h*k2/2);
k4=myfun(u1+h,v1+h*k3);
v2=v1+h*(k1+2*k2+2*k3+k4)/6;
u2=u1+h;
u1=u2;
v1=v2;
err=abs(y2-v2);
3、应用举例
4、测试结果
请输入初始点x1=1
请输入初始条件y1=2
请输入终止条件xn=6
请输入初始步长h1=0.1
h=0.10000000,x=1.00000000,y=2.00000000
h=0.02500000,x=1.02500000,y=2.02515952
h=0.02500000,x=1.05000000,y=2.05065134
1.3偏微分方程求解
参考附件:
《偏微分方程(PDEs)的MATLAB数值解法.pdf》。
1.4差分方程求解
《Z变换和差分方程的Matlab求解.pdf》。
1.5本章小结
本章主要介绍了代数方程、常微分方程、偏微分方程和差分方程的求解方法和示例讲解,这也是我在前几年培训过程中遇到的常用方法,希望大家牢记。
给大家推荐个轻音乐——夜的钢琴曲,石进的,晚上写教程的时候听的,感觉不错。
第二章微分和积分求解
本章主要介绍利用MATLAB自带的工具箱函数来实现微分、积分和差分的求解方法。
为了区别系统函数和自己编写的函数,我在有的函数后面加了_t以便区分。
2.1一般差分diff
使用diff可以对一组数据进行差分,如一组序列
,那么它的一阶差分为:
这就是差分的简单概念,具体XX相关资料参考。
调用方式:
Y=diff(X)
Y=diff(X,n)
Y=diff(X,n,dim)
其中,X是待差分序列,n差分的阶数(在微分方程里就是微分阶数),dim是所需要差分的纬度,Y是差分的结果。