数值积分微分方程Word文件下载.docx

数值积分微分方程Word文件下载.docx

《数值积分微分方程Word文件下载.docx》由会员分享，可在线阅读，更多相关《数值积分微分方程Word文件下载.docx（18页珍藏版）》请在冰豆网上搜索。

ifN==1

fprintf（'

Datahasonlyoneinterval'

）

return;

end

ifN==2

I=h/3*（f

（1）+4*f

（2）+f（3））;

ifN==3

I=3/8*h*（f

（1）+3*f

（2）+3*f（3）+f（4））;

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;

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:

m））;

例题求

。

解先编制

的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）

00.15640.30900.45400.58780.70710.8090

0.89100.95110.98771.00000.98770.95110.8910

0.80900.70710.58780.45400.30900.15640.0000

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（'

程序功能用复化梯形公式求定积分值

functionI=Trapezd（f_name,a,b,n）

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）]）;

补充例题求

y=sin（x）;

I=Trapezd（’sin_x’,0,pi,4）

I=1.8961

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中的一个时间点。

[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

输出

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类型

特点

说明

非刚性

一步算法；

4，5阶Runge-Kutta方程；

累计截断误差达（△x）3

大部分场合的首选算法

2，3阶Runge-Kutta方程；

使用于精度较低的情形

多步法；

Adams算法；

高低精度均可到10-3～10-6

计算时间比ode45短

适度刚性

采用梯形算法

适度刚性情形

刚性

Gear’s反向数值微分；

精度中等

若ode45失效时，可尝试使用

一步法；

2阶Rosebrock算法；

低精度

当精度较低时，计算时间比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’时，返回相应的事件记录

OutputFcn

odeplot、odephas2、odephas3、odeprint

若无输出参量，则solver将执行下面操作之一：

画出解向量中各元素随时间的变化；

画出解向量中前两个分量构成的相平面图；

画出解向量中前三个分量构成的三维相空间图；

随计算过程，显示解向量

OutputSel

正整数向量

[]

若不使用缺省设置，则OutputFcn所表现的是那些正整数指定的解向量中的分量的曲线或数据。

若为缺省值时，则缺省地按上面情形进行操作

Refine

正整数k>

1

k=1

若k>

1，则增加每个积分步中的数据点记录，使解曲线更加的光滑

Jacobian

若为‘on’时，返回相应的ode函数的Jacobi矩阵

Jpattern

为‘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）];

再在命令窗口中执行：

MU=7;

Y0=[1;

[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;

（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;

plot（x,y）

程序文件命名为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

-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（'

x0,y0,xend,M）

M为步长数

function[X,Y]=Rungkt4（fxy,x0,y0,xend,M）

h=（xend-x0）/M;

X=zeros（1,M+1）;

Y=zero