微分和差分方程.docx
《微分和差分方程.docx》由会员分享，可在线阅读，更多相关《微分和差分方程.docx（34页珍藏版）》请在冰豆网上搜索。
微分和差分方程
Part6——微分和差分方程
这次教程有点多，望大家耐心看哦！
第一章几类方程求解
1.1代数方程求解
下面主要介绍几类常用的方程求根方法。
（1）多项式——roots
调用示例：
r=roots（p）
其中，p为多项式的系数，从高次到低次顺序，r为对应的多个根。
如多项式：
，求解程序如下：
>>p=[1-6-72-27];
>>r=roots（p）
r=
12.1229
-5.7345
-0.3884
注意：
通过其它函数也可以多项式的根，但是不能求出所有的根，而roots可以求出多项式所有的根。
下面将会看到。
参考函数：
poly,residue
（2）一元函数——fzero
调用示例：
x=fzero（fun,x0）
x=fzero（fun,x0,options）
[x,fval]=fzero（...）
其中，fun为待求函数的名称，x0为初始值，x为解，fval为最优解对应的值（对于求根来说，理想情况下为0）。
如：
，求根程序如下：
f=@（x）x.^3-2*x-5;
z=fzero（f,2）
再看看
（1）中的示例，程序
f=@（x）x.^3-6*x.^2-72*x-27;
z1=fzero（f,10）
z2=fzero（f,-4）
z3=fzero（f,0）
很容易发现，对于不同的初始值，最后的根也不一样，最后的结果分别为：
12.1229、-5.7345、-0.3884，而这三个分别为
（1）中roots得到的解。
实际上fzero函数的实现方式是以某个初始点来迭代得到方程解的，因为初始值的选取对最终结果影响很大。
下面说到的牛顿迭代之类的算法也是类似情况，对初始值的选取很敏感，而智能算法对初始值的选取并不敏感。
注意：
f=@（x）...是匿名函数的表达方式，即该函数没有名字，但是通过@来指定变量，或者待求变量。
参考函数：
fminbnd,optimset,function_handle（@）,AnonymousFunctions
小技巧：
快捷注释与注释消除键，选取待注释语句，然后ctrl+R为注释，ctrl+T为消除注释。
（3）非线性方程组——fsolve
调用示例：
x=fsolve（fun,x0）
x=fsolve（fun,x0,options）
[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
0.5671
fval=
1.0e-006*
-0.4059
-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的符号解
solve（'a*x^2+b*x+c','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）
%%设置精度
x=newton（'myfun',[1;2],50,0.0000001）
其中newton函数如下：
functionx=newton（funname,x0,Maxgen,tol）
%作用:
通过牛顿迭代算法求解非线性方程组
%调用方式：
x=newton（'f_name',x0）
%x=newton（'f_name',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;
end
fprintf（'n=%3.0f,x=%12.5e,y=%12.5e,',n,x,y）
ifn>M
fprintf（'\n'）;
disp（'Warning:
iterationsexceedsthelimitation,probabledevergent'）;
end
注意：
通过调节最大迭代次数很容易发现，一般在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）;
结果如下：
x=
0.66701.3340
fval=
-8.2258
exitflag=
1
注意：
exitflag是算法终止标志，1代表正常终止，其它参考help。
这里的函数参数传递是通过@gatestfun匿名方式的，而上面我们编写的newton是通过'myfun'字符串方式传递的。
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'）;
x0=[1,1];
[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
disp（'没有为方程创建名为myfun.m的函数文件,请参照下例建立它'）;
disp（'functionz=myfun（x,y）'）;
disp（'z=y-2*x/y;'）;
disp（'并将该文件保存在work文件夹下'）;
end
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
fprintf（'第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y）;
plot（X,Y,'o'）
holdon
end
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、算法原理
单从每一步看，步长越小，截断误差就越小，但随着步长的缩小，在一定求解范围内所要完成的步数就增加了；步数的增加不但引起计算量的增大，而且可能导致舍入误差的严重积累。
因此同积分的数值计算一样，微分方程的数值解法也有个选择步长的问题。
在选择步长时，需要考虑两个问题：
1>怎样衡量和检验计算结果的精度？
2>如何依据所获得的精度处理步长？
考察经典的四阶龙格-库塔公式（2.2），从节点
出发，先以
为步长求出一个近似值
，由于公式的局部截断误差为
，故有：
（2.3）
然后将步长折半，即取
为步长从
跨两步到
，再求得一个近似值
，每跨一步的截断误差是
，因此有：
（2.4）
比较（2.3）式和（2.4）式我们们可以看到，步长折半后，误差大约减少到
。
即有：
由此得到下列的估算公式：
这样，可以通过检查步长，折半前后两次计算结果的偏差：
来判定所选的步长是否合适。
具体的说，将区分以下两种情况处理：
1、对于给定的精度
，如果
，反复将步长折半进行计算，直至
为止；这时取最终得到的
作为结果。
2、如果
，反复将步长加倍，直到
为止，这时再将步长折半一次，就得到所要的结果。
这种通过加倍或折半处理步长的方法称为变步长方法。
2、程序实现
functionmain
%
%测试主函数
%
ifexist（'myfun.m'）==0
disp（'没有为方程创建名为myfun.m的函数文件,请参照下例建立它'）;
disp（'functionz=myfun（x,y）'）;
disp（'z=y-2*x/y;'）;
disp（'并将该文件保存在work文件夹下'）;
end
ifexist（'self.m'）==0
disp（'把work文件夹里没有self.m文件'）;
end
eps=10^（-8）;
x1=input（'请输入初始点x1='）;
y1=input（'请输入初始条件y1='）;
xn=input（'请输入终止条件xn='）;
h1=input（'请输入初始步长h1='）;
h=h1;
fprintf（'h=%10.8f,x=%10.8f,y=%10.8f\n',h,x1,y1）;%输出初始条件
ifh>abs（x1-xn）
disp（'初始步长取得过大，超过了求解区间的长度'）;
x1=xn+1;%终止运算
end
whilex1<=xn
[u2,v2,h,err]=self（x1,y1,h）;
H=h;
half_err=err;
double_err=err;
%当误差过大时，反复缩小步长
whilehalf_err>eps
H=h;
[u2,v2,h,err]=self（x1,y1,H）;
half_err=err;
end
%当误差过小时，不断将步长增大
whiledouble_errH=2*H;
[u2,v2,h,err]=self（x1,y1,H）;
double_err=err;
end
%误差合适，最后进行调整运算
ifdouble_err>=eps
H=H/2;
[u2,v2,h,err]=self（x1,y1,H）;
end
%输出此点结果,为下一节点的计算提供初始值
fprintf（'h=%10.8f,x=%10.8f,y=%10.8f\n',H,u2,v2）;
x1=u2;
y1=v2;
h=h1;
plot（x1,y1,'o'）
holdon
end
其中，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;
end
err=abs（y2-v2）;
3、应用举例
functionz=myfun（x,y）
%%微分方程为：
dy/dx=y-2*x/y
z=y-2*x/y;
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是差分的结果。
应用举例：
x=
升级会员 