数值积分微分方程.docx
- 文档编号:4105433
- 上传时间:2022-11-27
- 格式:DOCX
- 页数:18
- 大小:96.60KB
数值积分微分方程.docx
《数值积分微分方程.docx》由会员分享,可在线阅读,更多相关《数值积分微分方程.docx(18页珍藏版)》请在冰豆网上搜索。
数值积分微分方程
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=zero
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 积分 微分方程