微分方程数值解.docx

微分方程数值解.docx

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

微分方程数值解

微分方程数值解

4.1

当常微分方程能解析求解时，可利用Matlab符号工具箱中的功能找到精确解.见下例

求解方程,,,.键入：

yyy，，,20

symsxy%定义符号变量diff_equ=‘D2y+2*Dy-y=0’;%D2y表示,,,Dy=y,yy=dsolve（diff_equ,‘x’）%定义x为自变量y=cl*exp（（2^（1/2）-1）*x+c2*exp（-（2^（1/2）+1）*x）

%表达式中含c1与c2，表示通解.

%初始条件为y（0）=0,,y（0）=1时，按如下方式调用y=dsolve（diff_equ,‘y（0）=0’,‘Dy（0）=1’,‘x’）y=1/4*2^（1/2）*exp（（2^（1/2）-1）*x）—

1/4*2^（1/2）*exp（-（2^（1/2）+1）*x）

%画出函数y=y（x）的图形

ezplot（y,[-2,2]）

图形具体形式请上机试之.

在方程无法获得解析解的情况下，可方便地获得数值解.下面的例子说明用Matlab求

数值解的方法及应注意的问题.

[例1]求解范德堡（vanderpol）方程

2dxdx2,,，,,

（1）0xx2dtdt

求解高阶方程，必须等价地变换为一阶微分方程组，对本例，通过定义两个新的变量，

实现这一变换

yxydxdt1,2/,,则令dydty1/2,

2dydtyyy2/（11）*21,,,,

编写求解程序分为两部分，第一部分为待求解的方程，存盘的文件名为,待求解方程

的函数名.m,，第二部分为求解主程序，本例中取名为main1.m.

首先编写待求解方程的文件.文件存盘名为“vdpol.m”.M,

functionyprime=vdpol（,）ty

yprime

（1）=y

（2）;

;（1

（1）^2）*

（2）

（1）,,yyy

mu=2yprime=[yprime

（1）;yprime

（2）];yprime

（2）=mu*说明函数yprime=vdpol中.定义为自变量，的形式取决于求解方程的阶数，本（,）tyyt

例中，,为解向量，为导数向量.yprime，y

（2）yyyy,[

（1）,

（2）],

（1）

（1）

（1）,y

yprime,,，函数返回vanderpol方程的导数列向量.因为所求结果为方程数值解，

（2）

（1）,y

所以各向量维数只有在主程序求解时定下精度后才能确定.

主程序定名为main1.m，你可用你所喜欢的其它名子，但vdpol.m除外.clearfunctions

%调试程序时，放置这一语句是必要的.它清除前边已编译的存在于内存中的废弃程序[]=ode23（‘vdpol’,[0,30],[1,0]）;ty,

y1=y（:

1）;%解曲线.

y2=y（:

2）;%解曲线的导数.

polt（‘__’）tyty,1,,2,

说明龙格_库塔的2阶与4阶改进型求解公式的实现，其指令分别为：

[]=ode23（‘f’,tsx,0,options）tx,

[]=ode45（‘f’,tsx,0,options）tx,

其中可由系统依据精度要求自动设定，亦可由使用者依据实际需要自己确定，分别说明之.ts

（1）若令tstttf,[0,1,,]，则输出在指定时刻tttf0,1,,给出，当tstktf,0:

:

时，

输出在区间[0,]ttf的等分点上给出，为步长.k

（2）若tsttft,[0,],0为自变量初值，tf为终值，此时，options决定自变量的维数，t中的时间点不是等间隔的，这是为了保证所需的相对精度，积分算法改变了步长.用于t

3,6设定误差限的参数options可缺省，此时系统设定相对误差为，绝对误差为，若1010

自行设定误差限，可用如下语句：

options=odeset（‘reltol’,,‘abstol’,）rtat

这里的与分别为设定的相对与绝对误差.rtat

须注意的是无论用哪种方法确定ttf0,的取值方式，必须由使用者确定且应与相匹配.x0t

y0[1,0],ts,[0,30]y（0）1,y（0）0,为初始条件，本例中，因为，这意味着解曲线，，x0

一般说，当解nnn个未知函数的方程组时，为维向量，共含有个初始条件.x0

两个输出参数是列向量xx与矩阵，它们具有相同的行数，而矩阵的列数等于方程t

组的个数，本例中y（:

1）y（:

2）的列数为2，其中，为自变量上各点函数值，为上各ytt

点导数值.

最后，提请读者注意的是：

ode45也不总是比ode23好，在很多时候，低阶算法更有

效，有关微分方程数值解法的更进一步信息，请参考数值分析方面的书籍.有些参考书提

供了一些关于算法选择和如何处理那些时间常数变化范围大的病态方程的非常实用的算法.4.2-

设有一阶方程与初始条件

yfxy,（,）,（4.1）,yxy（）,00,

其中适当光滑，关于满足Lipschitz条件，即存在使fxy（,）Ly

fxyfxyLyy（,）（,）,,,1212

则（4.1）式的解存在且惟一.

关于yyx,（）的解析解一般难以求到或根本无解析解，因而，实际问题中，通常，采

用差分的方法.在一系列离散点xxx,,,,yyy,,,,上寻求其数值近似解.12n12n

相邻两个节点间的间距xxnhn,，,,1,2,hxx,,称为步长，一般地取等步长，则hn0nnn，1

1、欧拉方法

在区间[,]xx上用差商nn，1

yxyx（）（）,nn，1h代替（4.1）式中,[,]xxxxxy，对fxy（,）中在上取值还是，而形成向前欧拉公式nn，1nn，1

与向后欧拉公式.

（1）向前欧拉公式

xfxy（,）取左端点，得如下公式n

yxyxhfxyx（）（）（,（））,，（4.2）nnnn，1

从yxy（）,x点出发，由初值代入（4.2）求得000

yyhfxy,，（,）（4.3）1000

反复利用（4.2），有

yyhfxyn,，,（,）0,1,2,（4.4）nnnn，1

其几何意义如图4.1所示.

y图中yyx,（）为方程（4.1）的精确PP43P2解曲线，其上任意点（,）xy处切线斜率为误差P1yyx,（）32Pxy（,）fxy（,）.从初值点出发，用该P0000y0点斜率fxy（,）xx,作一直线段，在001yyx,（）yx（）3处得到Pxy（,）y，由（4.2）式确定，1111y3再从Pxy（,）fxy（,）出发，以为斜率11111

作直线段，在xx,Pxy（,）处得到，2222xxxxxxO03124

PPP，012

作为积分曲线yx（）的近似，用图4.1yyx,（）n，1这一过程继续下去，形成折线表示在xy处的精确值，为解的近似值，不难得到n，1n，1

2h32,,（）（）（）（）yxyyxOhOh,,，,nnn，，112

P，1这一误差称为局部截断误差.若一种算法局部截断误差为Oh（），则称该算法具有阶P

精度，所以向前欧拉公式具有1阶精度.

（2）向后欧拉公式

若[,]xxxx中取中的，则有如下公式：

fxy（,）nn，1n，1

yyhfxyn,，,（,）0,1,2,（4.5）nnnn，，，111

称式（4.5）为向后欧拉公式，因为此式中y未知，故称其为隐式公式，无法用其直n，1接计算y，一般用向前欧拉公式产生初值.n，1

（0）yyhfxyn,，,（,）0,1,2,11nnnn，，

再按下式迭代

（1）（）kk，yyhfxykn,，,,（,）,0,1,,0,1,nnnn，，，111

其误差估计如下

2h32,,（）（）（）（）yxyyxohoh,,，,nnn，，112

精度亦为1阶，将向前欧拉公式（4.4）与向后欧拉公式（4.5）及它们的误差的几何说明作一对比，是十分有益的，见图4.2.

y为讨论局部截断误差，在图4.2中设点

APxy（,）落在积分曲线yyx,（）上，按式nnn

yyx,（）（4.4）及式（4.5）分别得,P点为与，ABn，1B且PAB,yyx,（）点一定在积分曲线上相应n点的上、下两边,所以将式（4.4）与（4.5）,

平均之，一定能得到更好的结果.

xxx（3）梯形公式nn，1将向前与向后欧拉公式加以平均得到所图4.2谓梯形公式

hyyfxyfxyn,，，,[（,）（,）]0,1,2,（4.6）nnnnnn，，，1112

3其局部截断误差为Oh（），具有2阶精度.

（4）改进的欧拉公式

为使计算简单，又免去迭代的繁复，将公式（4.6）简化为两步

yyhfxy,，（,）nnnn，1

h（4.7）yyfxyfxyn,，，,[（,）（,）],0,1,2,nnnnn，，11n，12或写为

h,yykk,，，（）nn，112,2,1nn（4.8）n,0,1,2,kfxy,（,）,

211nn，kfxyhk,，（,）,

最后指出，上述欧拉方法可推广至微分方程组，如

yfxyz,（,,）,

,zgxyz,（,,）,,yxy（）,00,

zxz（）,,00向前欧拉公式为

yyhfxyz,，（,,）,nnnnn，1n,0,1,2,,zzhgxyz,，（,,）,nnnnn，1

2、龙格_库塔方法

由微分中值定理

[（）（）]/（）,01yxyxhyxh,,，,,,,nnn，1

又因为,,yxhfxhyxh（）（,（）），,，，,,,yfxy,（,），所以nnn从而有

yxyxhfxhyxh（）（）（）,（）,，，，,,（4.9）nnnn，1令[,]xx，称其为区间上的平均斜率，由（4.9）可知，给kfxhyxh,，，（,（））,,nn，1nn

出一种平均斜率，可相应导出一种算法.向前欧拉公式中，精度低.改进欧kfxy,（,）nn

1拉公式中取[,]xxkfxyfxy,，（（,）（,）），精度提高，下面，我们在区间内nn，1nnn，1n，12

多取几个点，将其斜率加权平均，就能构造出精度更高的计算公式，公式的推导不再具体

给出，只开列具体结果.

（1）2阶龙格_库塔公式

yyhkk,，，（）,,,nn，11122,kfxy,（,）（4.10）1nn,

21nnkfxahyhka,，，,,（,）,0,1,,,

1,其中,,,，,,,1,,1a，由于4个未知数只有3个方程，所以解不惟一，若令1222a

1，即得改进的欧拉公式，具有2阶精度.,,,，,,,,1a122

（3）4阶龙格_库塔公式

只给出精典格式中最常用的一种.

h,yykkkk,，，，，（22）nn，11234,6,kfxy,（,）,1nn

hh,（4.11）kfxyk,，，（,）,21nn22,

hh,kfxyk,，，（,）32nn,22,kfxhyhk,，，（,）,43nn,

其计算精度为4阶

4.3

1、模型与问题

[例2]单摆运动

图4.3中一根长的细线，一端固定，另一端悬挂质量为l

m的小球，在重力作用下，小球处于竖直的平衡位置.现使小球偏离平衡位置一个小的角度，然后使其自由运动，在不,

考虑空气阻力情形下，小球将沿弧线作周期一定的简谐运动.

为平衡位置，在小球摆动过程中，当与平衡位置夹,,0

角为,mgsin,时，小球所