实验报告七 常微分方程初值问题的数值解法
常微分方程初值问题的数值解法
贵州师范大学数学与计算机科学学院学生实验报告课程名称： 数值分析 班级： 实验日期： 年 月 日 学 号： 姓名： 指导教师： 实验成绩： 一、实验名称实验六: 常微分方程初值问题数值解法 二、实验目的及要求1. 让学生掌握用Euler 法, Runge-Kutta 法求解常微分方程初值问题.2. 培养Matlab 编程与上机调试能力. 三、实验环境每人一台计算机，要求安装Windows XP 操作系统，Microsoft office2003、MATLAB6.5(或7.0). 四、实验内容1. 取步长h=0.1,0.05,0.01, ,用Euler 法及经典4阶Runge-Kutta 法求解初值问题⎩⎨⎧=≤≤++-=1)0()10(2222'y t t t y y 要求:1) 画出准确解(准确解22t e y t +=-)的曲线,近似解折线;2) 把节点0.1和0.5上的精确解与近似解比较,观察误差变化情况.2. 用 Euler 法,隐式Euler 法和经典4阶R-K 法取不同步长解初值问题⎪⎩⎪⎨⎧=∈-=21)0(],1,0[,50'y x y y 并画出曲线观察稳定性. 注：题1必须写实验报告五、算法描述及实验步骤Euler 法：输入 000),(,,,),,(y a x x h b a y x f = 输出 Euler 解y步1 ),,2,1(;m n h n a x hab m n =⨯+=-⇐步2 对1,,2,1,0-=m n 执行),(1n n n n y x f h y y ⨯+⇐+ 步3 输出T m y y y y ),,,(21 = 经典4阶R-K 法：输入 000),(,,,),,(y a x x h b a y x f = 输出 4阶R-K 解y 步1 ),,2,1(;m n h n a x hab m n =⨯+=-⇐步2 对1,,2,1,0-=m n 执行),(1n n y x f K ⇐，)5.0,(15.02hK y x f K n n +⇐+， )5.0,(25.03hK y x f K n n +⇐+，),(314hK y x f K n n +⇐+ )22(643211K K K K hy y n n ++++⇐+ 步3 输出T m y y y y ),,,(21 =六、调试过程及实验结果>> shiyan6Y1 =0.8000 0.6620 0.5776 0.5401 0.5441 0.5853 0.6602 0.7662 0.9009 1.0627 Y2 =0.8287 0.7103 0.6388 0.6093 0.6179 0.6612 0.7366 0.8419 0.9753 1.1353e1 =0.0287e2 =4.2469e-006 e1 =0.0738e2 =1.1609e-005注：至于h=0.05、0.01的情况将程序中的h值作相应的改动即可得。
常微分方程初值问题的数值解法
1 1 2 1 , 2 p 2
这里有 3 个未知 数， 2 个方程。
存在无穷多个解。所有满足上式的格式统称为2阶龙格 - 库 塔格式。注意到，p 1, 1 2 1 就是改进的欧拉法。
Step 1: 将 K2 在 ( xi , yi ) 点作 Taylor 展开
K 2 f ( xi ph, yi phK1 ) f ( xi , yi ) phf x ( xi , yi ) phK1 f y ( xi , yi ) O( h2 )
y( xi ) phy( xi ) O(h2 )
d f ( x, y) dx 首先希望能确定系数 1、2、p，使得到的算法格式有 2阶 dy 精度，即在 yi y( xi ) 的前提假设下，使得 f x ( x, y) f y ( x, y) dx Ri y( xi 1 ) yi 1 O(h3 ) f x ( x, y) f y ( x, y) f ( x, y) y( x )
y( x0 ) y0 yk 1 yk h f ( xk , yk 1 )
, k 0,1,...
隐式欧拉法的求解: 利用迭代的思路进行.
yi 1 yi hf ( xi , yi 1 )
变换为
y
( k 1) i 1
yi hf ( xi , y )
y i 1 K1 K2
1 1 y i h K 1 K 2 2 2 f ( xi , yi ) f ( xi h, yi hK 1 )
步长一定是一个h 吗？
§2 Runge-Kutta Method
常微分方程初值问题数值解法
常微分方程初值问题的数值解法在自然科学、工程技术、经济和医学等领域中，常常会遇到一阶常微分方程初值问题：(,),,(),y f x y a x b y a y '=≤≤⎧⎨=⎩ (1) 此处f 为,x y 的已知函数，0y 是给定的初始值。
本章讨论该问题的数值解法，要求f 在区域{(,)|,}G x y a x b y =≤≤<∞内连续，并对y 满足Lipschitz 条件，从而初值问题(1)有唯一的连续可微解()y y x =，且它是适定的。
1 几个简单的数值积分法1.1 Euler 方法（1）向前Euler 公式（显式Euler 公式）10(,),0,1,2,,(),n n n n y y hf x y n y y a +=+=⎧⎨=⎩(2) 其中h 为步长。
由此便可由初值0y 逐步算出一阶常微分方程初值问题(1)的解()y y x =在节点12,,x x 处的近似值12,,y y 。
该公式的局部截断误差为2()O h ，是一阶方法。
（2）向后Euler 公式（隐式Euler 公式）1110(,),0,1,2,,(),n n n n y y hf x y n y y a +++=+=⎧⎨=⎩(3) 这是一个隐格式，也是一阶方法。
这类隐格式的计算比显格式困难，一般采用迭代法求解。
首先用向前Euler 公式提供迭代初值，然后迭代计算：(0)1(1)()111(,),(,),0,1,2,n n n n k k n n n n y y hf x y y y hf x y k +++++⎧=+⎨=+=⎩ (4)1.2 梯形方法1110[(,)(,)],2(),(0,1,2,)n n n n n n h y y f x y f x y y y a n +++⎧=++⎪⎨⎪=⎩= (5) 这也是一个隐格式，是二阶方法。
一般也采用迭代法求解。
迭代公式如下：(0)1(1)()111(,),[(,)(,)],0,1,2,2n n n n k k n n n n n n y y hf x y h y y f x y f x y k +++++⎧=+⎪⎨=++=⎪⎩ (6)1.3 改进的Euler 方法11110(,),[(,)(,)],0,1,2,,2(),n n n n n n n n n n y y hf x y h y y f x y f x y n y y a ++++⎧=+⎪⎪=++=⎨⎪=⎪⎩(7) 为了便于上机编程计算，(7)可改写为110(,),(,),0,1,2,,1(),2(),p n n n cn n p n p c y y hf x y y y hf x y n y y y y y a ++=+⎧⎪=+⎪⎪=⎨=+⎪⎪=⎪⎩(8) 该格式是显式，也是二阶方法。
常微分方程的数值解法实验报告
常微分方程的数值解法专业班级：信息软件 姓名：吴中原 学号：120108010002 一、实验目的1、熟悉各种初值问题的算法，编出算法程序；2、明确各种算法的精度与所选步长有密切关系；通过计算更加了解各种 算法的优越性。
二、实验题目1、根据初值问题数值算法，分别选择二个初值问题编程计算；2、试分别取不同步长，考察某节点j x处数值解的误差变化情况； 3、试用不同算法求解某初值问题，结果有何异常； 4、分析各个算法的优缺点。
三、实验原理与理论基础（一） 欧拉法算法设计对常微分方程初始问题(6-1)(6-2)用数值方法求解时，我们总是认为(6-1)、(6-2)的解存在且唯一。
欧拉法是解初值问题的最简单的数值方法。
从(6-2)式由于y (x 0) = y 0已给定，因而可以算出),()('000y x f x y =。
设x 1 = h 充分小，则近似地有：),()(')()(00001y x f x y hx y x y =≈-(6-3)记 ,n ,,i x y y i i 10 )(== 从而我们可以取),(0001y x hf y y ==作为)(1x y 的近似值。
利用1y 及f (x 1, y 1)又可以算出)(2x y 的近似值：),(1112y x hf y y +=一般地，在任意点()h n x n 11+=+处)(x y 的近似值由下式给出),(1n n n n y x hf y y +=+(6-4)这就是欧拉法的计算公式，h 称为步长。
⎪⎩⎪⎨⎧==)( ),(d d 00y x y y x f x y（二）四阶龙格-库塔法算法设计：欧拉公式可以改写为：()111,i i i i y y k k hf x y +=+⎧⎪⎨=⎪⎩，它每一步计算(),f x y 的值一次，截断误差为()2o h 。
改进的欧拉公式可以改写为：()()()11212112,,i i i i i i y y k k k hf x y k hf x h y k +⎧=++⎪⎪=⎨⎪=++⎪⎩，它每一步要计算(),f x y 的值两次，截断误差为()3o h 。
第7章 常微分方程初值问题的数值解法
例1 函数f ( t , y ) = t y 在区域D0 = {( t , y ) | 1 ≤ t ≤ 2, −3 ≤ y ≤ 4}
关于y满足Lipschitz条件，相应的Lipschitz常数可取为L = 2
3 存在性定理 定理1 设函数f ( t , y )在凸集D ⊂ R 2中有定义，若存在常数
(7.2.7)
称为显式Runge-Kutta(龙格-库塔 )方法，简称R-K方法，
其中正整数N 称为R-K方法的级，所有ci , ai , bij 都是待定 常数。
根据定义(7.2.7)，N 级R-K方法(7.9)的局部截断误差为
Rn+1 = y( t n+1 ) − y( t n ) − h∑ ci ki
dy 其斜率为 = f ( t0 , y0 ) dt ( t0 , y0 ) 由 点 斜 式 写 出 切线 方 程 dy y = y0 + ( t − t0 ) = y0 + ( t − t0 ) f ( t0 , y0 ) dt ( x0 , y )
0
等步长为h,则t1 - t0 = h, 可由切线算出 y1 ： 则 y1 = y0 + hf ( t0 , y0 ) 按此逐步计算y( tn ), 在tn +1处的值 ： yn+1 = yn + hf ( tn , yn ) y 注意： 这是“ 注意 ： 这是 “ 折 yN 线法” 而非“ 线法 ” 而非 “ 切 线法” 线法 ” 除第一个 点是曲线切线外， 点是曲线切线外 ， 其他点不是切线 y2 而是折线(如右 y1 y0 图所示)。 图所示 。
பைடு நூலகம்
则称数值解法(7.5)为显式方法。否则，称数值解法(7.3) 为隐式方法。
常微分方程初值问题
实验报告专用纸实验项目名称常微分方程初值问题的数值解课程名称计算机数值方法教师评语及成绩：实验成绩：教师签字：（请按照实验报告的有关要求书写，一般必须包括：1、实验目的；2、实验内容；3、实验步骤与方法；4、实验数据与程序清单；5、出现的问题及解决方法；6、实验结果、结果分析与体会等内容。
）1、实验目的（1）在计算机上实现Euler法及Runge-Kutta法求解一阶常微分方程初值问题；（2）通过求常微分方程初值问题的数值解，加深对其数值方法及有关基础理论知识的理解。
2、实验内容（1）用Euler法求解常微分方程初值问题；（2）用四阶Runge-Kutta求解常微分方程初值问题。
3、实验原理（1）Euler法：（2）改进的Euler法：（3）Runge-Kutta法：4、实验步骤（1）复习相关知识，为编写程序做准备；（2）编写调式程序（3）用调式好的程序解决实际问题（4）记录相关数据（5）进行结果分析（6）根据要求撰写论文5、实验程序%Ordinary_differential.mformat long;%步长h=1/20;%右边界xf=2;%设置步长x=0:h:xf;%已知函数Func=@(x,y)(y-((2*x)/y));%Euler%初始化y数组y_E=zeros(1,length(x));y_E(1)=1;%复现公式&迭代计算for i=1:(length(x)-1)y_E(i+1)=y_E(i)+Func(x(i),y_E(i))*h; end%Improve Euler%初始化y数组y_IE=zeros(1,length(x));y_IE(1)=1;%复现公式&迭代计算for i=1:(length(x)-1)p=y_IE(i)+h*Func(x(i),y_IE(i));c=y_IE(i)+h*Func(x(i+1),p);y_IE(i+1)=(1/2)*(p+c);end%Rung Kutta%初始化y数组y_RK=zeros(1,length(x));y_RK(1)=1;%复现公式&迭代计算for i=1:(length(x)-1)k1=Func(x(i),y_RK(i));k2=Func(x(i)+0.5*h,y_RK(i)+0.5*h*k1);k3=Func((x(i)+0.5*h),(y_RK(i)+0.5*h*k2));k4=Func((x(i)+h),(y_RK(i)+k3*h));y_RK(i+1)=y_RK(i)+(1/6)*(k1+2*k2+2*k3+k4) *h;end%Exact%使用ode45计算精确函数x0=0;y0=1;yspan=[x0xf];[x_ode45,y_ode45]=ode45(Func,yspan,y0);%绘图subplot(221)plot(x_ode45,y_ode45,'-');xlabel('x');ylabel('y');legend('Exact');subplot(222)plot(x,y_E,'-');xlabel('x');ylabel('y');legend('Euler');subplot(223)plot(x,y_IE,'-');xlabel('x');ylabel('y');legend('Improve Euler');subplot(224)plot(x,y_RK,'-');xlabel('x');ylabel('y');legend('Rung Kutta');%tableres=[round(transpose(x),8),round(transpose(y_E),8), round(transpose(y_IE),8),round(transpose(y_RK),8),round(y_ode45,8)];table=array2table(res,'VariableNames',{'x','Euler','Improve Euler','Rung Kutta','Exact'});table;6、实验结果及分析（1）实验结果table=41×5tablex Euler Improve Euler Rung Kutta Exact_______________________________________________ 011110.05 1.05 1.04886905 1.04880886 1.048810180.1 1.0977381 1.09556112 1.09544514 1.095445380.15 1.14351536 1.14034459 1.14017546 1.140174580.2 1.18757368 1.18343698 1.183216 1.183216060.25 1.23011131 1.22501749 1.22474493 1.22474526 0.3 1.27129351 1.26523585 1.26491113 1.26491111 0.35 1.31126017 1.30421877 1.30384056 1.30384029 0.4 1.3501313 1.34207454 1.34164088 1.34164095 0.45 1.38801111 1.3788967 1.37840498 1.37840509 0.5 1.42499118 1.41476663 1.41421368 1.41421362 0.55 1.4611528 1.44975574 1.44913781 1.44913769 0.6 1.49656893 1.48392711 1.48323985 1.48323991 0.65 1.53130567 1.51733679 1.51657526 1.51657529 0.7 1.56542352 1.55003492 1.54919353 1.54919345 0.75 1.59897836 1.58206653 1.58113904 1.58113897 0.8 1.63202233 1.61347233 1.61245178 1.61245182 0.85 1.66460451 1.64428925 1.64316793 1.64316793 0.9 1.69677156 1.67455096 1.67332034 1.673320260.95 1.72856823 1.7042883 1.70293895 1.702938891 1.76003786 1.73352962 1.73205115 1.73205117 1.05 1.79122279 1.76230109 1.76068206 1.76068204 1.1 1.82216476 1.79062697 1.78885479 1.78885471 1.15 1.85290524 1.81852981 1.81659066 1.8165906 1.2 1.8834858 1.84603071 1.84390938 1.84390939 1.25 1.91394844 1.87314942 1.87082923 1.87082919 1.3 1.94433584 1.89990456 1.89736719 1.89736709 1.35 1.97469176 1.92631373 1.92353905 1.92353898 1.4 2.00506125 1.95239363 1.94935958 1.94935956 1.45 2.03549101 1.9781602 1.97484254 1.97484247 1.5 2.06602967 2.0036287 2.00000085 2.00000073 1.55 2.09672813 2.0288138 2.0248466 2.0248465 1.6 2.12763984 2.05372973 2.04939117 2.04939111 1.65 2.15882113 2.07839027 2.07364525 2.073645151.72.19033159 2.10280892 2.09761891 2.097618771.752.22223435 2.12699893 2.12132168 2.121321541.82.2545965 2.15097339 2.14476252 2.144762411.852.28748942 2.17474529 2.16794994 2.167949781.92.3209892 2.19832764 2.19089198 2.190891781.952.35517701 2.22173348 2.21359628 2.213596082 2.39013954 2.24497601 2.23607008 2.2360699（2）结果分析根据数据及图像可以看出。
7 常微分方程初值问题的数值解法
h
2
y
( p)
(xn )
h
p
p! (8.7)
y
( p1)
( n )
( p 1 )!
其中
y , f ( x, y) ,, ' ' , ' ' y f x ( x, y) f x ( x, y) y f x f y f ( p) y
,
y ( n )
,
2! y ( n )
,,
( x n1 x n ) 或 (8 .2 )
2
y ( x n ) h f ( x n , y ( x n )) yn1 yn h f ( x n , yn )
h
2
2
( n 0 , 1, 2 , )
(8 .2 ) 叫 E u ler 显 格 式 ， 可 循 环 求 解 。 方 法 二 、 数 值 微 分 法 — — 实 质 还 是 T a ylo r 展 开 法 （ 略 ）
p1
).
定理 8 . 1
设 f ( x , y ) 关于第 2 个变元满足 Lipschitz 条件 | f ( x 1 , y 1 ) f ( x 2 , y 2 ) | L | y 1 , y 2 | 且 Euler 显格式的局部截断误差 | n | e
(ba )
满足 | R n | 1 ]h
' '
[ ]
(8 .9 )
p! 假 设 y n＝ y ( x n )， 则 （ 8 .7） 和 （ 8 .9） 得 局 部 截 断 误 差 为 Rn1 y ( x n1 ) yn1 h
p1
方法三 将方程
实验报告七常微分方程初值问题的数值解法
实验报告七常微分方程初值问题的数值解法浙江大学城市学院实验报告课程名称数值计算方法实验项目名称常微分方程初值问题得数值解法实验成绩指导老师(签名) 日期2０15/１2/16一、实验目得与要求1．用Ｍａｔlab软件掌握求微分方程数值解得欧拉方法与龙格-库塔方法;2．通过实例学习用微分方程模型解决简化得实际问题。
二、实验内容与原理编程题2－1要求写出Matlab源程序（m文件),并有适当得注释语句;分析应用题2－２,２-3,2-4,2－5要求将问题得分析过程、Ｍatlab源程序与运行结果与结果得解释、算法得分析写在实验报告上。
2-1编程编写用向前欧拉公式与改进欧拉公式求微分方程数值解得Ｍａtlab程序,问题如下:在区间内个等距点处，逼近下列初值问题得解,并对程序得每一句添上注释语句。
Euler法 y＝euleｒ(ａ,b，ｎ,y0，f,f1，ｂ1)改进Euler法y=eulerpro（ａ,b,ｎ,ｙ0，f,f1,b１)2-2分析应用题假设等分区间数,用欧拉法与改进欧拉法在区间内求解初值问题并作出解得曲线图形,同时将方程得解析解也画在同一张图上,并作比较,分析这两种方法得精度。
2-3分析应用题用以下三种不同得方法求下述微分方程得数值解,取画出解得图形,与精确值比较并进行分析。
1)欧拉法；２)改进欧拉法;3)龙格－库塔方法;2-4分析应用题考虑一个涉及到社会上与众不同得人得繁衍问题模型。
假设在时刻(单位为年）,社会上有人口人,又假设所有与众不同得人与别得与众不同得人结婚后所生后代也就是与众不同得人。
而固定比例为得所有其她得后代也就是与众不同得人。
如果对所有人来说出生率假定为常数,又如果普通得人与与众不同得人得婚配就是任意得,则此问题可以用微分方程表示为:其中变量表示在时刻社会上与众不同得人得比例,表示在时刻人口中与众不同得人得数量。
1）假定与,当步长为年时,求从到解得近似值,并作出近似解得曲线图形。
2)精确求出微分方程得解,并将您当时在分题(b)中得到得结果与此时得精确值进行比较。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
课程名称 数值计算方法实验项目名称 常微分方程初值问题的数值解法实验成绩 指导老师（签名 ） 日期 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)10y t y t y '=-⎧⎨=⎩ 并作出解的曲线图形，同时将方程的解析解也画在同一张图上，并作比较，分析这两种方法的精度。
2-3 分析应用题用以下三种不同的方法求下述微分方程的数值解，取10h =201(0)1y 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 =114，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 分别为设定的相对误差和绝对误差。
常用选项见下表。
选项名 功能 可选值 省缺值 AbsTol 设定绝对误差 正数 16e - RelTol 设定相对误差 正数 13e - InitialStep 设定初始步长 正数 自动 MaxStep 设定步长上界正数 10tspan MaxOrder 设定ode15s 的最高阶数 1,2,3,4,5 5 Stats 显示计算成本统计on,off off BDF设定ode15s 是否用反向差分on,offoff例：解微分方程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 yans =……>> plot(t ,y ,‘o-’,) %解函数图形表示>> 45(,[0,4],1)ode odefun %不用输出变量，则直接输出图形>> [],45(,0:4,1)t y ode odefun =；[],t yans =三. 操作方法与实验步骤（包括实验数据记录和处理）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));endplot(x,y)hold on% 求微分方程的精确解x1=linspace(a,b,100);'精确解为's=dsolve(f1,b1,'x')syms xy1=zeros(1,100);fori=1:100y1(i)=subs(s,x,x1(i)); endplot(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;fori=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); endplot(x,y)hold on% 求微分方程的精确解x1=linspace(a,b,100); '精确解为's=dsolve(f1,b1,'x') syms xy1=zeros(1,100); for i=1:100y1(i)=subs(s,x,x1(i)); endplot(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 8Columns 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 80Columns 81 through 88Columns 89 through 96Columns 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改进欧拉法的精度比向前欧拉法更高。