实验报告七 常微分方程初值问题的数值解法

课程名称 数值计算方法

实验项目名称 常微分方程初值问题的数值解法

实验成绩 指导老师（签名 ） 日期 2015/12/16

一. 实验目的和要求

1． 用Matlab 软件掌握求微分方程数值解的欧拉方法和龙格－库塔方法； 2． 通过实例学习用微分方程模型解决简化的实际问题。

二. 实验内容和原理

编程题2-1要求写出Matlab 源程序(m 文件)，并有适当的注释语句；分析应用题2-2，2-3，2-4，2-5要求将问题的分析过程、Matlab 源程序和运行结果和结果的解释、算法的分析写在实验报告上。

2-1 编程

编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab 程序，问题如下：

在区间[],a b 内(1)N +个等距点处，逼近下列初值问题的解，并对程序的每一句添上注释语句。

0(,)()y f x y a x b y a y '=≤≤=

Euler 法 y=euler(a,b,n,y0,f,f1,b1)

改进Euler 法 y=eulerpro(a,b,n,y0,f,f1,b1)

2-2 分析应用题

假设等分区间数100n =，用欧拉法和改进欧拉法在区间[0,10]t ∈内求解初值问题

()()20

(0)10

y t y t y '=-??

=? 并作出解的曲线图形，同时将方程的解析解也画在同一张图上，并作比较，分析这两种方法的精度。

2-3 分析应用题

用以下三种不同的方法求下述微分方程的数值解，取10h =

201

(0)1

y y x x y '=+≤≤??

=?

画出解的图形，与精确值比较并进行分析。 1）欧拉法； 2）改进欧拉法； 3）龙格－库塔方法；

2-4 分析应用题

考虑一个涉及到社会上与众不同的人的繁衍问题模型。假设在时刻t (单位为年)，社会上有人口()x t 人，又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人。而固定比例为r 的所有其他的后代也是与众不同的人。如果对所有人来说出生率假定为常数b ，又如果普通的人和与众不同的人的婚配是任意的，则此问题可以用微分方程表示为：

()

(1())dp t rb p t dt

=- 其中变量()()()i p t x t x t =表示在时刻t 社会上与众不同的人的比例，()i x t 表示在时刻t 人口中与众不同的人的数量。

1）假定(0)0.01,0.02p b ==和0.1r =，当步长为1h =年时，求从0t =到50t =解()p t 的近似值，并作出近似解的曲线图形。

2）精确求出微分方程的解()p t ，并将你当50t =时在分题(b)中得到的结果与此时的精确值进行比较。

【MATLAB 相关函数】

求微分方程的解析解及其数值的代入

dsolve(‘egn1’， ‘egn2’，L ‘x ’) subs (expr, {x,y,…}, {x1,y1,…} )

其中‘egn i ’表示第i 个方程，‘x ’表示微分方程中的自变量，默认时自变量为t 。 subs 命令中的expr 、x 、y 为符合型表达式，x 、y 分别用数值x1、x2代入。 >> syms x y z

>> subs('x+y+z',{x,y,z},{1,2,3})

ans =

6 >> syms x

>> subs('x^2',x,2)

ans =

4

>> s=dsolve(‘12Dy y ∧=+’, ‘(0)1y =’, ‘x ’) ans =

tan(14)x pi -*

>> syms x

>> subs(s,x,2)

ans =

右端函数(,)f x y 的自动生成

f= inline(‘expr ’, ’var1’, ‘var2’,……)

其中’expr ’表示函数的表达式，’var1’, ‘var2’ 表示函数表达式中的变量，运

行该函数，生成一个新的函数表达式为f (var1, var2, ……)。 >> f=inline('x+3*y','x','y') f =

Inline function: f(x,y) = x+3*y >> f(2,3)

ans =

11

4，5阶龙格－库塔方法求解微分方程数值解

[t,x]=ode45(f,ts,x0,options)

其中f 是由待解方程写成的m 文件名；x0为函数的初值；t,x 分别为输出的自变量和函

数值(列向量)，t 的步长是程序根据误差限自动选定的。若ts=[t0,t1,t2,…,tf]，则输出在自变量指定值，等步长时用ts=t0:k:tf ，输出在等分点；options 用于设定误差

限(可以缺省，缺省时设定为相对误差3

10-，绝对误差6

10-)，程序为：

options=odeset(‘reltol ’,rt,’abstol ’,at)，这里rt,at 分别为设定的相对误差和绝对误差。常用选项见下表。

选项名 功能 可选值 省缺值 AbsTol 设定绝对误差 正数 16e - RelTol 设定相对误差 正数 13e - InitialStep 设定初始步长 正数 自动 MaxStep 设定步长上界

正数 10tspan MaxOrder 设定ode15s 的最高阶数 1,2,3,4,5 5 Stats 显示计算成本统计

on,off off BDF

设定ode15s 是否用反向差分

on,off

off

例：解微分方程

204(0)1t y y t y y ?'

=-<

??=?

在命令窗口执行

>> odefun = inline (‘2*y t y -’, ‘t ’, ‘y ’)；

>> [],45(,[0,4],1)t y ode odefun =； >> [],t y

ans =

……

>> plot(t ,y ,‘o-’,) %解函数图形表示

>> 45(,[0,4],1)ode odefun %不用输出变量，则直接输出图形

>> [],45(,0:4,1)t y ode odefun =；[],t y

ans =

三. 操作方法与实验步骤（包括实验数据记录和处理）

2-1 编程

编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab 程序，问题如下：

在区间[],a b 内(1)N +个等距点处，逼近下列初值问题的解，并对程序的每一句添上注释语句。

0(,)()y f x y a x b y a y '=≤≤=

Euler 法 y=euler(a,b,n,y0,f,f1,b1)

改进Euler 法 y=eulerpro(a,b,n,y0,f,f1,b1)

Euler 法

y=euler(a,b,n,y0,f,f1,b1) y=zeros(1,n+1);

y(1)=y0;

h=(b-a)/n;

x=a:h:b;

for i=1:n;

y(i+1)=y(i)+h*f(x(i),y(i));

end

plot(x,y)

hold on

% 求微分方程的精确解

x1=linspace(a,b,100);

'精确解为'

s=dsolve(f1,b1,'x')

syms x

y1=zeros(1,100);

for

i=1:100

y1(i)=subs(s,x,x1(i)); end

plot(x1,y1,'r')

title('红色代表精确解')

改进Euler法

y=eulerpro(a,b,n,y0,f,f1,b1) % 求微分方程的数值解

y=zeros(1,n+1);

y(1)=y0;

h=(b-a)/n;

x=a:h:b;

for

i=1:n;

T1=f(x(i),y(i));

T2=f(x(i+1),y(i)+h*T1);

y(i+1)=y(i)+(h/2)*(T1+T2); end

plot(x,y)

hold on

% 求微分方程的精确解

x1=linspace(a,b,100); '精确解为'

s=dsolve(f1,b1,'x') syms x

y1=zeros(1,100); for i=1:100

y1(i)=subs(s,x,x1(i)); end

plot(x1,y1,'r')

title('红色代表精确解')

2-2分析应用题

假设等分区间数100n =，用欧拉法和改进欧拉法在区间[0,10]t ∈内求解初值问题

()()20

(0)10y t y t y '=-??

=?

并作出解的曲线图形，同时将方程的解析解也画在同一张图上，并作比较，分析这两种方法

的精度。

(1)向前欧拉法

>> euler(0,10,100,10,inline('y-20','x','y'),'Dy=y-20','y(0)=10') ans =

精确解为 s =

20 - 10*exp(x) ans =

+005 *

Columns 1 through 8

Columns 9 through 16 Columns 17 through 24 Columns 25 through 32 Columns 33 through 40 Columns 41 through 48 Columns 49 through 56 Columns 57 through 64 Columns 65 through 72 Columns 73 through 80

Columns 81 through 88

Columns 89 through 96

Columns 97 through 101

（2）改进欧拉法

>> eulerpro(0,10,100,10,inline('y-20','x','y'),'Dy=y-20','y(0)=10') ans =

精确解为

s =

20 - 10*exp(x)

ans =

+005 *

Columns 1 through 8 Columns 9 through 16 Columns 17 through 24 Columns 25 through 32 Columns 33 through 40 Columns 41 through 48 Columns 49 through 56 Columns 57 through 64 Columns 65 through 72 Columns 73 through 80 Columns 81 through 88 Columns 89 through 96 Columns 97 through 101

改进欧拉法的精度比向前欧拉法更高。

2-3分析应用题

用以下三种不同的方法求下述微分方程的数值解，取10h = 201

(0)1

y y x x y '=+≤≤??

=?

画出解的图形，与精确值比较并进行分析。 1）欧拉法； 2）改进欧拉法；

2-4分析应用题

考虑一个涉及到社会上与众不同的人的繁衍问题模型。假设在时刻t (单位为年)，社会上有人口()x t 人，又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人。而固定比例为r 的所有其他的后代也是与众不同的人。如果对所有人来说出生率假定为常数b ，又如果普通的人和与众不同的人的婚配是任意的，则此问题可以用微分方程表示为：

()

(1())dp t rb p t dt

=- 其中变量()()()i p t x t x t =表示在时刻t 社会上与众不同的人的比例，()i x t 表示在时刻t 人口中与众不同的人的数量。

1）假定(0)0.01,0.02p b ==和0.1r =，当步长为1h =年时，求从0t =到50t =解()p t 的近似值，并作出近似解的曲线图形。

2）精确求出微分方程的解()p t ，并将你当50t =时在分题(b)中得到的结果与此时的精确值

进行比较。

1）

>>

euler(0,50,50,,inline('','t','p'),'Dp=','p(0)= 1') ans =

精确解为

s =

1 - 99/(100*exp(x/500)) ans =

Columns 1 through 8 Columns 9 through 16 Columns 17 through 24 Columns 25 through 32 Columns 33 through 40 Columns 41 through 48 Columns 49 through 51

（2）

>> dsolve('Dp=','p(0)=','t') ans =1 - 99/(100*exp(t/500))

>> 1 - 99/(100*exp)

ans =

与欧拉法求得的精确值差0,0001四. 实验结果与分析