浙江大学城市学院实验报告
课程名称 数值计算方法
实验项目名称 常微分方程初值问题的数值解法
实验成绩 指导老师（签名 ） 日期 2015/12/16
一. 实验目的和要求
1． 用Matlab软件掌握求微分方程数值解的欧拉方法和龙格－库塔方法； 2． 通过实例学习用微分方程模型解决简化的实际问题。
二. 实验内容和原理
编程题2-1要求写出Matlab源程序(m文件)，并有适当的注释语句；分析应用题2-2，2-3，2-4，2-5要求将问题的分析过程、Matlab源程序和运行结果和结果的解释、算法的分析写在实验报告上。
2-1 编程
编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab程序，问题如下： 在区间?a,b?内(N?1)个等距点处，逼近下列初值问题的解，并对程序的每一句添上注释语句。
y??f(x,y) a?x?b y(a)?y0
Euler法 y=euler(a,b,n,y0,f,f1,b1)
改进Euler法 y=eulerpro(a,b,n,y0,f,f1,b1)
2-2 分析应用题
假设等分区间数n?100，用欧拉法和改进欧拉法在区间t?[0,10]内求解初值问题
?y?(t)?y(t)?20 ??y(0)?10并作出解的曲线图形，同时将方程的解析解也画在同一张图上，并作比较，分析这两种方法的精度。
2-3 分析应用题
用以下三种不同的方法求下述微分方程的数值解，取h?10 ??y??y?2x?y(0)?10?x?1
画出解的图形，与精确值比较并进行分析。 1）欧拉法； 2）改进欧拉法； 3）龙格－库塔方法；
2-4 分析应用题
考虑一个涉及到社会上与众不同的人的繁衍问题模型。假设在时刻t(单位为年)，社会上有人口x(t)人，又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人。而固定比例为r的所有其他的后代也是与众不同的人。如果对所有人来说出生率假定为常数b，又如果普通的人和与众不同的人的婚配是任意的，则此问题可以用微分方程表示为：
dp(t)?rb(1?p(t)) dt其中变量p(t)?xi(t)x(t)表示在时刻t社会上与众不同的人的比例，xi(t)表示在时刻t人口
中与众不同的人的数量。
1）假定p(0)?0.01,b?0.02和r?0.1，当步长为h?1年时，求从t?0到t?50解p(t)的近似值，并作出近似解的曲线图形。
2）精确求出微分方程的解p(t)，并将你当t?50时在分题(b)中得到的结果与此时的精确值进行比较。
【MATLAB相关函数】
? 求微分方程的解析解及其数值的代入
dsolve(‘egn1’， ‘egn2’，? ‘x’) subs (expr, {x,y,…}, {x1,y1,…} )
其中‘egni’表示第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(‘Dy?1?y?2’, ‘y(0)?1’, ‘x’) ans =
tanx(?1?4pi )>> syms x >> subs(s,x,2)
ans =
-0.3721
? 右端函数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用于设定误差
限(可以缺省，缺省时设定为相对误差10，绝对误差10)，程序为：
options=odeset(‘reltol’,rt,’abstol’,at)，这里rt,at分别为设定的相对误差和绝对误差。常用选项见下表。
选项名 AbsTol RelTol InitialStep MaxStep MaxOrder Stats BDF
例：解微分方程
功能 设定绝对误差 设定相对误差 设定初始步长 设定步长上界 设定ode15s的最高阶数 显示计算成本统计 设定ode15s是否用反向差分 可选值 正数 正数 正数 正数 1,2,3,4,5 on,off on,off 省缺值 1e?6 1e?3 自动 tspan10 5 off off ?3?6t??0?t?4?y?y?2 y ?
?y(0)?1?在命令窗口执行
?? odefun = inline(‘y?2*ty’, ‘t’, ‘y’)；
?? ?t,y??ode45(odefun,[0,4],1)； ?? ?t,y?
ans =
0 1.0000 0.0502 1.0490 0.1005 1.0959 0.1507 1.1408 ……
3.8507 2.9503 3.9005 2.9672 3.9502 2.9839 4.0000 3.0006
?? plot(t,y,‘o-’,) %解函数图形表示
?? ode45(odefun,[0,4],1) %不用输出变量，则直接输出图形
?? ?t,y??ode45(odefun,0:4,1)；?t,y?
ans =
0 1.0000 1.0000 1.7321 2.0000 2.2361 3.0000 2.6458 4.0000 3.0006
三. 操作方法与实验步骤（包括实验数据记录和处理）
2-1 编程
编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab程序，问题如下： 在区间?a,b?内(N?1)个等距点处，逼近下列初值问题的解，并对程序的每一句添上注释语句。
y??f(x,y) a?x?b y(a)?y0
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);