数值积分与微分方程.docx
《数值积分与微分方程.docx》由会员分享，可在线阅读，更多相关《数值积分与微分方程.docx（19页珍藏版）》请在冰点文库上搜索。
数值积分与微分方程
2.3数值积分
2.3.1一元函数的数值积分
函数1quad、quadl、quad8
功能数值定积分，自适应Simpleson积分法。
格式q=quad（fun,a,b）%近似地从a到b计算函数fun的数值积分，误差为10-6。
若给fun输入向量x，应返回向量y，即fun是一单值函数。
q=quad（fun,a,b,tol）%用指定的绝对误差tol代替缺省误差。
tol越大，函数计算的次数越少，速度越快，但结果精度变小。
q=quad（fun,a,b,tol,trace,p1,p2,…）%将可选参数p1,p2,…等传递给函数fun（x,p1,p2,…），再作数值积分。
若tol=[]或trace=[]，则用缺省值进行计算。
[q,n]=quad（fun,a,b,…）%同时返回函数计算的次数n
…=quadl（fun,a,b,…）%用高精度进行计算，效率可能比quad更好。
…=quad8（fun,a,b,…）%该命令是将废弃的命令，用quadl代替。
例2-40
>>fun=inline（‘3*x.^2./（x.^3-2*x.^2+3）’）;equivalentto:
functiony=funn（x）
y=3*x.^2./（x.^3-2*x.^2+3）;
>>Q1=quad（fun,0,2）
>>Q2=quadl（fun,0,2）
计算结果为：
Q1=
3.7224
Q2=
3.7224
补充：
复化simpson积分法程序
程序名称Simpson.m
调用格式I=Simpson（'f_name',a,b,n）
程序功能用复化Simpson公式求定积分值
输入变量f_name为用户自己编写给定函数
的M函数而命名的程序文件名
a为积分下限
b为积分上限
n为积分区间
划分成小区间的等份数
输出变量I为定积分值
程序
functionI=simpson（f_name,a,b,n）
h=（b-a）/n;
x=a+（0:
n）*h;
f=feval（f_name,x）;
N=length（f）-1;
ifN==1
fprintf（'Datahasonlyoneinterval'）
return;
end
ifN==2
I=h/3*（f
（1）+4*f
（2）+f（3））;
return;
end
ifN==3
I=3/8*h*（f
（1）+3*f
（2）+3*f（3）+f（4））;
return;
end
I=0;
if2*floor（N/2）==N
I=h/3*（2*f（N-2）+2*f（N-1）+4*f（N）+f（N+1））;
m=N-3;
else
m=N;
end
I=I+（h/3）*（f
（1）+4*sum（f（2:
2:
m））+2*f（m+1））;
ifm>2
I=I+（h/3）*2*sum（f（3:
2:
m））;
end
例题求
。
解先编制
的M函数。
程序文件命名为sin_x.m。
functiony=sin_x（x）
y=sin（x）
将区间4等份，调用格式为
I=Simpson（’sin_x’,0,pi,4）
计算结果为
y=
00.70711.00000.70710.0000
I=
2.0046
将区间20等份，调用格式为
I=Simpson（’sin_x’,0,pi,20）
计算结果为
y=
00.15640.30900.45400.58780.70710.8090
0.89100.95110.98771.00000.98770.95110.8910
0.80900.70710.58780.45400.30900.15640.0000
I=
2.0000
重做上例2—40：
simpson（'funn',0,2,100）
函数2trapz
功能梯形法数值积分
格式T=trapz（Y）%用等距梯形法近似计算Y的积分。
若Y是一向量，则trapz（Y）为Y的积分；若Y是一矩阵，则trapz（Y）为Y的每一列的积分；若Y是一多维阵列，则trapz（Y）沿着Y的第一个非单元集的方向进行计算。
T=trapz（X,Y）%用梯形法计算Y在X点上的积分。
若X为一列向量，Y为矩阵，且size（Y,1）=length（X），则trapz（X,Y）通过Y的第一个非单元集方向进行计算。
T=trapz（…,dim）%沿着dim指定的方向对Y进行积分。
若参量中包含X，则应有length（X）=size（Y,dim）。
例2-41
>>X=-1:
.1:
1;
>>Y=1./（1+25*X.^2）;
>>T=trapz（X,Y）
计算结果为：
T=
0.5492
补充：
复化梯形积分法程序
程序名称Trapezd.m
调用格式I=Trapezd（'f_name',a,b,n）
程序功能用复化梯形公式求定积分值
输入变量f_name为用户自己编写给定函数
的M函数而命名的程序文件名
a为积分下限
b为积分上限
n为积分区间
划分成小区间的等份数
输出变量I为定积分值
程序
functionI=Trapezd（f_name,a,b,n）
h=（b-a）/n;
x=a+（0:
n）*h;
f=feval（f_name,x）;
I=h*（sum（f）-（f
（1）+f（length（f）））/2）;
hc=（b-a）/100;
xc=a+（0:
100）*hc;
fc=feval（f_name,xc）;
plot（xc,fc,'r'）;
holdon;
title（'TrapezoidalRule'）;xlabel（'x'）;ylabel（'y'）;
plot（x,f）;
plot（x,zeros（size（x）））;
fori=1:
n;
plot（[x（i）,x（i）],[0,f（i）]）;
end
补充例题求
。
解先编制
的M函数。
程序文件命名为sin_x.m。
functiony=sin_x（x）
y=sin（x）;
将区间4等份，调用格式为
I=Trapezd（’sin_x’,0,pi,4）
计算结果为
I=1.8961
将区间20等份，调用格式为
I=Trapezd（’sin_x’,0,pi,20）
计算结果为
I=1.9959
图A.5表示了复化梯形求积的过程。
（1）区间4等份
（2）区间20等份
重做上例2-41：
functiony=li2_41（x）
y=1./（1+25*x.^2）;
I=Trapezd（’li2_41’,-1,1,100）
函数3rat,rats
功能有理分式近似。
虽然所有的浮点数值都是有理数，有时用简单的有理数字（分子与分母都是较小的整数）近似地表示它们是有必要的。
函数rat将试图做到这一点。
对于有连续出现的小数的数值，将会用有理式近似表示它们。
函数rats调用函数rat，且返回字符串。
格式[N,D]=rat（X）%对于缺省的误差1.e-6*norm（X（:
）,1），返回阵列N与D，使N./D近似为X。
[N,D]=rat（X,tol）%在指定的误差tol范围内，返回阵列N与D，使N./D近似为X。
rat（X）、rat（X…）%在没有输出参量时，简单地显示x的连续分数。
S=rats（X,strlen）%返回一包含简单形式的、X中每一元素的有理近似字符串S，若对于分配的空间中不能显示某一元素，则用星号表示。
该元素与X中其他元素进行比较而言较小，但并非是可以忽略。
参量strlen为函数rats中返回的字符串元素的长度。
缺省值为strlen=13，这允许在78个空格中有6个元素。
S=rats（X）%返回与用MATLAB命令formatrat显示 X相同的结果给S。
例2-42
>>s=1-1/2+1/3-1/4+1/5-1/6+1/7
>>formatrat
>>S1=rats（s）
>>S2=rat（s）
>>[n,d]=rat（s）
>>PI1=rats（pi）
>>PI2=rat（pi）
计算结果为：
s=
0.7595
S1=
319/420
S2=
1+1/（-4+1/（-6+1/（-3+1/（-5））））
n=
319
d=
420
PI1=
355/113
PI2=
3+1/（7+1/（16））
2.3.2二元函数重积分的数值计算
函数1dblquad
功能矩形区域上的二重积分的数值计算
格式q=dblquad（fun,xmin,xmax,ymin,ymax）%调用函数quad在区域[xmin,xmax,ymin,ymax]上计算二元函数z=f（x,y）的二重积分。
输入向量x，标量y，则f（x,y）必须返回一用于积分的向量。
q=dblquad（fun,xmin,xmax,ymin,ymax,tol）%用指定的精度tol代替缺省精度10-6，再进行计算。
q=dblquad（fun,xmin,xmax,ymin,ymax,tol,method）%用指定的算法method代替缺省算法quad。
method的取值有@quadl或用户指定的、与命令quad与quadl有相同调用次序的函数句柄。
q=dblquad（fun,xmin,xmax,ymin,ymax,tol,method,p1,p2,…）%将可选参数p1,p2,..等传递给函数fun（x,y,p1,p2,…）。
若tol=[]，method=[]，则使用缺省精度和算法quad。
例2-43
>>fun=inline（’y./sin（x）+x.*exp（y）’）;
>>Q=dblquad（fun,1,3,5,7）
计算结果为：
Q=
3.8319e+003
2.4常微分方程数值解
函数ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb
功能常微分方程（ODE）组初值问题的数值解
参数说明：
solver为命令ode45、ode23,ode113,ode15s,ode23s,ode23t,ode23tb之一。
Odefun为显式常微分方程y’=f（t,y），或为包含一混合矩阵的方程M（t,y）*y’=f（t,y）。
命令ode23只能求解常数混合矩阵的问题；命令ode23t与ode15s可以求解奇异矩阵的问题。
Tspan积分区间（即求解区间）的向量tspan=[t0,tf]。
要获得问题在其他指定时间点t0,t1,t2,…上的解，则令tspan=[t0,t1,t2,…,tf]（要求是单调的）。
Y0包含初始条件的向量。
Options用命令odeset设置的可选积分参数。
P1,p2,…传递给函数odefun的可选参数。
格式[T,Y]=solver（odefun,tspan,y0）%在区间tspan=[t0,tf]上，从t0到tf，用初始条件y0求解显式微分方程y’=f（t,y）。
对于标量t与列向量y，函数f=odefun（t,y）必须返回一f（t,y）的列向量f。
解矩阵Y中的每一行对应于返回的时间列向量T中的一个时间点。
要获得问题在其他指定时间点t0,t1,t2,…上的解，则令tspan=[t0,t1,t2,…,tf]（要求是单调的）。
[T,Y]=solver（odefun,tspan,y0,options）%用参数options（用命令odeset生成）设置的属性（代替了缺省的积分参数），再进行操作。
常用的属性包括相对误差值RelTol（缺省值为1e-3）与绝对误差向量AbsTol（缺省值为每一元素为1e-6）。
[T,Y]=solver（odefun,tspan,y0,options,p1,p2…）将参数p1,p2,p3,..等传递给函数odefun，再进行计算。
若没有参数设置，则令options=[]。
1．求解具体ODE的基本过程：
（1）根据问题所属学科中的规律、定律、公式，用微分方程与初始条件进行描述。
F（y,y’,y’’,…,y（n）,t）=0
y（0）=y0,y’（0）=y1,…,y（n-1）（0）=yn-1
而y=[y;y
（1）;y
（2）;…,y（m-1）]，n与m可以不等
（2）运用数学中的变量替换：
yn=y（n-1）,yn-1=y（n-2）,…,y2=y1=y，把高阶（大于2阶）的方程（组）写成一阶微分方程组：
，
（3）根据
（1）与
（2）的结果，编写能计算导数的M-函数文件odefile。
（4）将文件odefile与初始条件传递给求解器Solver中的一个，运行后就可得到ODE的、在指定时间区间上的解列向量y（其中包含y及不同阶的导数）。
2．求解器Solver与方程组的关系表见表2-3。
表2-3
函数指令
含义
函数
含义
求解器
Solver
ode23
普通2-3阶法解ODE
odefile
包含ODE的文件
ode23s
低阶法解刚性ODE
选项
odeset
创建、更改Solver选项
ode23t
解适度刚性ODE
odeget
读取Solver的设置值
ode23tb
低阶法解刚性ODE
输出
odeplot
ODE的时间序列图
ode45
普通4-5阶法解ODE
odephas2
ODE的二维相平面图
ode15s
变阶法解刚性ODE
odephas3
ODE的三维相平面图
ode113
普通变阶法解ODE
odeprint
在命令窗口输出结果
3．因为没有一种算法可以有效地解决所有的ODE问题，为此，MATLAB提供了多种求解器Solver，对于不同的ODE问题，采用不同的Solver。
表2-4不同求解器Solver的特点
求解器Solver
ODE类型
特点
说明
ode45
非刚性
一步算法；4，5阶Runge-Kutta方程；累计截断误差达（△x）3
大部分场合的首选算法
ode23
非刚性
一步算法；2，3阶Runge-Kutta方程；累计截断误差达（△x）3
使用于精度较低的情形
ode113
非刚性
多步法；Adams算法；高低精度均可到10-3～10-6
计算时间比ode45短
ode23t
适度刚性
采用梯形算法
适度刚性情形
ode15s
刚性
多步法；Gear’s反向数值微分；精度中等
若ode45失效时，可尝试使用
ode23s
刚性
一步法；2阶Rosebrock算法；低精度
当精度较低时，计算时间比ode15s短
ode23tb
刚性
梯形算法；低精度
当精度较低时，计算时间比ode15s短
4．在计算过程中，用户可以对求解指令solver中的具体执行参数进行设置（如绝对误差、相对误差、步长等）。
表2-5Solver中options的属性
属性名
取值
含义
AbsTol
有效值：
正实数或向量
缺省值：
1e-6
绝对误差对应于解向量中的所有元素；向量则分别对应于解向量中的每一分量
RelTol
有效值：
正实数
缺省值：
1e-3
相对误差对应于解向量中的所有元素。
在每步（第k步）计算过程中，误差估计为：
e（k）<=max（RelTol*abs（y（k））,AbsTol（k））
NormControl
有效值：
on、off
缺省值：
off
为‘on’时，控制解向量范数的相对误差，使每步计算中，满足：
norm（e）<=max（RelTol*norm（y）,AbsTol）
Events
有效值：
on、off
为‘on’时，返回相应的事件记录
OutputFcn
有效值：
odeplot、odephas2、odephas3、odeprint
缺省值：
odeplot
若无输出参量，则solver将执行下面操作之一：
画出解向量中各元素随时间的变化；
画出解向量中前两个分量构成的相平面图；
画出解向量中前三个分量构成的三维相空间图；
随计算过程，显示解向量
OutputSel
有效值：
正整数向量
缺省值：
[]
若不使用缺省设置，则OutputFcn所表现的是那些正整数指定的解向量中的分量的曲线或数据。
若为缺省值时，则缺省地按上面情形进行操作
Refine
有效值：
正整数k>1
缺省值：
k=1
若k>1，则增加每个积分步中的数据点记录，使解曲线更加的光滑
Jacobian
有效值：
on、off
缺省值：
off
若为‘on’时，返回相应的ode函数的Jacobi矩阵
Jpattern
有效值：
on、off
缺省值：
off
为‘on’时，返回相应的ode函数的稀疏Jacobi矩阵
Mass
有效值：
none、M、
M（t）、M（t,y）
缺省值：
none
M：
不随时间变化的常数矩阵
M（t）：
随时间变化的矩阵
M（t,y）：
随时间、地点变化的矩阵
MaxStep
有效值：
正实数
缺省值：
tspans/10
最大积分步长
注意：
（1）求微分方程数值解的函数命令中，函数odefun必须以dx/dt为输出量，以t,x为输入量。
（2）用于解n个未知函数的方程组时，M函数文件中的待解方程组应以x的向量形式写成。
例题A.7解微分方程
，其中
首先，将导数表达式的右端编写成一个liA7.m函数程序：
functionyy=liA7（x,y）
yy=sin（x）;
然后直接调用：
[x,y]=ode23（'liA7',[0,pi],-1）
plot（x,y）
例4.用微分方程的数值解法求解微分方程
，设自变量t的初始值为0,终值为3pi,初始条件y（0）=0,y’（0）=0
解：
将高阶微分方程化为一阶微分方程组，即用变量代换：
=
这样，将导数表达式的右端编写成一个exf.m函数程序
functionxdot=exf（t,x）
u=1-（t.^2）/（2*pi）;
xdot=[01;-10]*x+[01]'*u;
然后，在主程序中调用已有的数值积分函数进行积分：
clf;t0=0;tf=3*pi;x0=[0;0]
[t,x]=ode23（'exf',[t0,tf],x0）
y=x（:
1）
例5，求二阶微分方程
的数值解
解：
首先变量代换：
这样，将导数表达式的右端编写成一个jie3.m函数程序
functionf=jie3（x,z）
f=[01;1/（2*x^2）-1-1/x]*z;
然后，在主程序中调用已有的数值积分函数进行积分：
[x,z]=ode23（'jie3',[pi/2,pi],[2;-2/pi]）
plot（x,z（:
1））
例2-45求解描述振荡器的经典的VerderPol微分方程
y（0）=1，y’（0）=0
令x1=y，x2=dy/dt，则
dx1/dt=x2
dx2/dt=μ（1-（x1）^2）*x2-x1
编写函数文件verderpol.m：
functionxprime=verderpol（t,x）
globalMU
xprime=[x
（2）;MU*（1-x
（1）^2）*x
（2）-x
（1）];
再在命令窗口中执行：
>>globalMU
>>MU=7;
>>Y0=[1;0]
>>[t,x]=ode45（‘verderpol’,0,40,Y0）;
>>x1=x（:
1）;x2=x（:
2）;
>>plot（t,x1,t,x2）
图形结果为图2-20。
图2-20VerderPol微分方程图
A.4.1改进的Euler法程序
程序名称Eulerpro.m
调用格式[X,Y]=Eulerpro（'fxy',x0,y0,xend,h）
程序功能解常微分方程
输入变量fxy为用户编写给定函数
的M函数文件名
x0,xend为起点和终点
y0为已知初始值
h为步长
输出变量X为离散的自变量
Y为离散的函数值
程序
function[x,y]=Eulerpro（fxy,x0,y0,xend,h）
n=fix（（xend-x0）/h）;
y
（1）=y0;
x
（1）=x0;
fork=2:
n
x（k）=0;
y（k）=0;
end
fori=1:
（n-1）
x（i+1）=x0+i*h;
y1=y（i）+h*feval（fxy,x（i）,y（i））;
y2=y（i）+h*feval（fxy,x（i+1）,y1）;
y（i+1）=（y1+y2）/2;
end
plot（x,y）
例题A.7解微分方程
，其中
。
解先编制
的M函数。
程序文件命名为fxy.m。
functionZ=fxy（x,y）
Z=sin（x）;
取步长0.1，调用格式为
[X,Y]=Eulerpro（‘fxy’,0,-1,pi,0.1）
计算结果如图A.6所示。
图A.6微分方程求解结果
x=
00.10000.20000.30000.40000.50000.60000.7000
0.80000.90001.00001.10001.20001.30001.40001.5000
1.60001.70001.80001.90002.00002.10002.20002.3000
2.40002.50002.60002.70002.80002.90003.0000
y=
-1.0000-0.9950-0.9801-0.9554-0.9211-0.8777-0.8255-0.7650
-0.6970-0.6219-0.5407-0.4541-0.3629-0.2681-0.1707-0.0715
0.02830.12790.22620.32220.41500.50360.58720.6649
0.73590.79960.85530.90250.94060.96930.9883
A.4.2Runge-Kutta法程序
程序名称RungKt4.m
调用格式[X,Y]＝RungKt4（'fxy',x0,y0,xend,M）
程序功能解常微分方程
输入变量fxy为用户编写给定函数
的M函数文件名
x0,xend为起点和终点
y0为已知初始值
M为步长数
输出变量X为离散的自变量
Y为离散的函数值
程序
function[X,Y]=Rungkt4（fxy,x0,y0,xend,M）
h=（xend-x0）/M;
X=zeros（1,M+1）;
Y=zer
升级会员 