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

实验报告七常微分方程

初值问题的数值解法 Document number【SA80SAB-SAA9SYT-SAATC-SA6UT-SA18】

浙江大学城市学院实验报告

课程名称 数值计算方法

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

实验成绩 指导老师（签名 ） 日期 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 +个等距点处，逼近下列初值问题的解，并对程序的每一句添上注释语句。

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)10y t y t y '=-⎧⎨=⎩

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

2-3 分析应用题

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

画出解的图形，与精确值比较并进行分析。

1）欧拉法；

2）改进欧拉法；

3）龙格－库塔方法；

2-4 分析应用题

考虑一个涉及到社会上与众不同的人的繁衍问题模型。假设在时刻t (单位为年)，

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

其中变量()()()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’， ‘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 =

>> 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 用于设定误差限(可以缺省，缺省时设定为相对误差310-，绝对误

差610-)，程序为：options=odeset(‘reltol ’,rt,’abstol ’,at)，这里rt,at

例：在命令窗口执行

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

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

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 +个等距点处，逼近下列初值问题的解，并对程序的每一句添上注释语句。

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