试验四常微分方程初值问题数值解法讲解

1、实 验 报 告课程名称 数值分析实验项目 常微分方程问题初值问题数值解法. 实验目的1、理解如何在计算机上实现用 Euler 法、改进 Euler 法、Runge Kutta 算法求一 阶常微分方程初值问题y （x） f （x, y）,x a,b y（a） y1的数值解。2、利用图形直观分析近似解和准确解之间的误差。3、学会 Matlab 提供的 ode45 函数求解微分方程初值问题。二、实验要求（ 1） 按照题目要求完成实验内容；（ 2） 写出相应的 Matlab 程序；（ 3） 给出实验结果 （可以用表格展示实验结果 ）；（ 4） 分析和讨论实验结果并提出可能的优化实验。（ 5） 写出实验

2、报告。三、实验步骤1、用编好的 Euler法、改进 Euler法计算书本 P167 的例 1、P171例题 3（1）取 h 0.1，求解初值问题y （x） y 2x,x 0,1,yy（0） 1（2）取 h 0.1，求解初值问题y （x） y x 1,x 0,0.5,y（0） 12、用 RungeKutta算法计算 P178例题、 P285实验任务（ 2）1)取 h 0.1，求解初值问题y (x) y2,x 0,0.5,y(0) 1(2)求初值问题y (x) 1 ( y x2 4x 1), x 0,0.5,2y(0) 0x 的解 y(x) 在 xi ih(h 0.05) 处的近似值 yi ，并与

3、问题的解析解 y(x) e 2 x2 1 相比较。3、用 Matlab 绘图函数 plot(x,y) 绘制 P285实验任务( 2)的精确解和近似解的图 形。4、使用 matlab中的 ode45求解 P285实验任务( 2)，并绘图。四、实验结果1、Euler 算法程序、改进 Euler 算法程序；Euler算法程序： function x,y=euler_f(ydot_fun, x0, y0, h, N) % Euler (向前)公式，其中%ydot_fun-%x0, y0 -%h -%N -%x -Xn%y -Yn一阶微分方程的函数初始条件区间步长区间的个数构成的向量构成的向量x=zer

4、os(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0; for n=1:Nx(n+1)=x(n)+h;y(n+1)=y(n)+h*feval(ydot_fun, x(n), y(n); end改进 Euler 算法程序：function x,y=euler_r(ydot_fun, x0, y0, h, N) 改进 Euler 公式，其中 ydot_fun - x0, y0 -h%一阶微分方程的函数初始条件区间步长 区间的个数 Xn 构成的向量% y - Yn 构成的向量x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=

5、y0; for n=1:Nx(n+1)=x(n)+h; ybar=y(n)+h*feval(ydot_fun, x(n), y(n); y(n+1)=y(n)+h/2*(feval(ydot_fun, x(n), y(n)+feval(ydot_fun, x(n+1), ybar); end2、用 Euler算法程序、改进 Euler算法求解 P167例题 1的运行结果；（1.）Euler 算法程序：Columns 1 through 80 0.1000 0.2000 0.3000 0.4000 0.50000.70000.6000Columns 9 through 110.8000 0.90

6、00 1.00001.50900.50001.4164Columns 1 through 81.0000 1.1000 1.1918 1.2774 1.3582 1.4351 1.5803Columns 9 through 111.6498 1.7178 1.7848（2）改进 Euler 算法：x =Columns 1 through 70 0.1000 0.2000 0.3000 0.4000 0.6000Columns 8 through 110.7000 0.8000 0.9000 1.0000y =Columns 1 through 71.0000 1.0959 1.1841 1.2

7、662 1.3434 1.4860Columns 8 through 111.5525 1.6165 1.6782 1.73793、RungeKutta 算法程序；function x,y=Runge_Kutta4(ydot_fun, x0, y0, h, N) % 标准四阶 Runge_Kutta 公式，其中 %ydot_fun - 一阶微分方程的函数%x0, y0 - 初始条件%h- 区间步长%N- 区间的个数% x - Xn 构成的向量% y - Yn 构成的向量x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0; for n=1:Nx(n+

8、1)=x(n)+h;k1=h*feval(ydot_fun, x(n), y(n);k2=h*feval(ydot_fun, x(n)+1/2*h, y(n)+1/2*k1);k3=h*feval(ydot_fun, x(n)+1/2*h, y(n)+1/2*k2);k4=h*feval(ydot_fun, x(n)+h, y(n)+k3); y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4);end4、用 RungeKutta 算法求解 P178例题、P285实验任务（2），计算结果如下（其 中 yi表示数值解， y（xi） 表示解析解，结果保留八位有效数字） ：求解 P17

9、8 例题：0 0.1000 0.2000 0.3000 0.4000 0.50001.0000 1.1111 1.2500 1.4286 1.6667 2.0000xi0.050.10.150.20.25yi-0.0222-0.0388-0.0498-0.0552-0.0550y(xi)-0.0222-0.0388-0.0498-0.0552-0.0550xi0.30.350.40.450.5yi-0.0493-0.0380-0.02130.00100.0288y(xi)-0.0493-0.0380-0.02130.00100.02885、P285 实验任务（ 2）精确解与近似解的图形比较6、

10、用 matlab 中的 ode45求解 P285 实验任务（ 2） ans =Columns 1 through 700.00010.00020.00030.00040.00090.00140-0.0001-0.0001-0.0002-0.0002-0.0005-0.0007Columns 8 through 140.00190.00240.00490.00740.00990.01250.0250-0.0010-0.0012-0.0024-0.0037-0.0049-0.0061-0.0118Columns 15 through 210.0375 0.0500 0.0625 0.0750 0.

11、0875 0.1000 0.1125-0.0172 -0.0222 -0.0268 -0.0312 -0.0351 -0.0388Columns 22 through 280.1250 0.1375 0.1500 0.1625 0.1750 0.1875-0.0450 -0.0475 -0.0497 -0.0516 -0.0532 -0.0543Columns 29 through 350.2125 0.2250 0.2375 0.2500 0.2625 0.2750-0.0556 -0.0558 -0.0556 -0.0550 -0.0541 -0.0528Columns 36 throug

12、h 420.3000 0.3125 0.3250 0.3375 0.3500 0.3625-0.0493 -0.0470 -0.0444 -0.0414 -0.0381 -0.0344Columns 43 through 490.3875 0.4000 0.4125 0.4250 0.4375 0.4500-0.0260 -0.0213 -0.0162 -0.0108 -0.0051 0.0010Columns 50 through 530.4718 0.4812 0.49060.50000.0125 0.0177 0.02320.0288-0.04200.2000-0.05520.2875-

13、0.05120.3750-0.03040.46250.0074五、讨论分析 误差一开始变大，然后变小，最后又变大六、改进实验建议可以通过提高保留的位数，使 Runge Kutta 算法结果更精确，也可以使数值解 和解析解的比较更加精确实验四 常微分方程初值问题数值解法这里只提供解答格式请同学自己按照上机实验报告格式填写）1、饿 uler 法求解初值问题function x,y=euler_f(ydot_fun, x0, y0, h, N) % Euler (向前)公式，其中% ydot_fun -% x0, y0 -%h-%N-%x- Xn%y- Yn一阶微分方程的函数初始条件区间步长区间的个

14、数构成的向量构成的向量x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0;for n=1:Nx(n+1)=x(n)+h;y(n+1)=y(n)+h*feval(ydot_fun, x(n), y(n);end用欧拉法计算 p167例 1ydot_fun=inline(y-2*x./y,x,y); x,y=euler_f(ydot_fun,0,1,0.1,10)x = Columns 1 through 110 0.100000000000000 0.200000000000000 0.300000000000000 0.50000000000000

15、00.600000000000000 0.700000000000000 0.800000000000000 1.000000000000000y =.000000000000000 1.1000000000000001.358212599560289 1.4351329186577961.508966253566332 1.5803382376552171.1918181818181821.6497834310477111.7847708324979820.4000000000000000.9000000000000001.2774378337147221.7177793478600872、

16、改进 Euler 公式function x,y=euler_r(ydot_fun, x0, y0, h, N)% 改进 Euler 公式，其中%ydot_fun - 一阶微分方程的函数%x0, y0 - 初始条件%h-区间步长%N-区间的个数%x-Xn 构成的向量%y-Yn 构成的向量x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0; for n=1:Nx(n+1)=x(n)+h; ybar=y(n)+h*feval(ydot_fun, x(n), y(n); y(n+1)=y(n)+h/2*(feval(ydot_fun, x(n), y(n)

17、+feval(ydot_fun, x(n+1), ybar); end用改进欧拉公式计算 P167 例题 1ydot_fun=inline(y-2*x./y,x,y); x,y=euler_r(ydot_fun,0,1,0.1,10)y =Columns 1 through 61.0000000000000001.343360151483999 Columns 7 through 11 1.485955602415669 1.7378674010354141.0959090909090911.4164019285369091.5525140913261461.1840965692429971.

18、6164747827520581.2662013608757761.678166363675186x =Columns 1 through 60 0.100000000000000 0.200000000000000 0.300000000000000 0.4000000000000000.500000000000000Columns 7 through 110.600000000000000 0.700000000000000 0.800000000000000 0.9000000000000001.000000000000000用欧拉法、改进欧拉法计算 P171 例题 3ydot_fun=

19、inline(-y+x+1,x,y); x,y=euler_f(ydot_fun,0,1,0.1,5)x =00.1000000000000000.2000000000000000.3000000000000000.400000000000000 0.500000000000000y =1.000000000000000 1.000000000000000 1.010000000000000 1.056100000000000 1.090490000000000 用改进欧拉法计算 P171例题 3 ydot_fun=inline(-y+x+1,x,y); x,y=euler_r(ydot_fu

20、n,0,1,0.1,5)x=0 0.100000000000000 0.200000000000000 0.400000000000000 0.5000000000000001.0290000000000000.300000000000000y=1.000000000000000 1.005000000000000 1.0190250000000001.070801950625000 1.1070757653156251.0412176250000003、标准四阶 Runge_Kutta 公式function x,y=Runge_Kutta4(ydot_fun, x0, y0, h, N) %

21、 标准四阶 Runge_Kutta 公式，其中 %ydot_fun - 一阶微分方程的函数%x0, y0 - 初始条件% h - 区间步长%N- 区间的个数%x- Xn 构成的向量% y - Yn 构成的向量x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0; for n=1:Nx(n+1)=x(n)+h;k1=h*feval(ydot_fun, x(n), y(n);k2=h*feval(ydot_fun, x(n)+1/2*h, y(n)+1/2*k1);k3=h*feval(ydot_fun, x(n)+1/2*h, y(n)+1/2*k2)

22、; k4=h*feval(ydot_fun, x(n)+h, y(n)+k3); y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4); end用四阶龙格库塔计算 P178 例题 ydot_fun=inline(y2,x,y); x,y=Runge_Kutta4(ydot_fun,0,1,0.1,5)x =0 0.100000000000000 0.200000000000000 0.400000000000000 0.5000000000000000.300000000000000y =1.000000000000000 1.111110490052194 1.24999799

23、2047015 1.4285661863014451.666653257250323 1.999963258950669用龙格库塔方法计算 P285实验任务( 2)ydot_fun=inline(-y+x2+4*x-1)./2,x,y); x,y=Runge_Kutta4(ydot_fun,0,0,0.05,10)x =Columns 1 through 6 0 0.050000000000000 0.250000000000000 Columns 7 through 11 0.300000000000000 0.500000000000000y = Columns 1 through 6 0

24、 -0.022190087076823 -0.055003093175772Columns 7 through 110.1000000000000000.1500000000000000.2000000000000000.3500000000000000.4000000000000000.450000000000000-0.038770573733692-0.049756511058554-0.055162578526671-0.049292018554665 -0.038042973450910-0.0212692404030080.0010162259975860.028800791009

25、418xx=0:0.05:0.5; z=exp(-xx./2)+xx.2-1; hold on plot(x,y,o);plot(xx,z,*)hold off 其图形如图一所示0.030.020.010-0.01-0.02-0.03-0.04-0.050 0.1 0.2 0.3 0.4 0.5-0.06图一 精确解与龙格库塔计算结果比较5、 用 Matlab 提供的 ode45求解 P285 实验任务( 2)ydot_fun=inline(-y+x2+4*x-1)./2,x,y); x,y=ode45(ydot_fun, 0,0.5, 0); x;yplot(x,y, * )ans = Co

26、lumns 1 through 60 0.000100475457260 0.0002009509145210.0003014263717810.0004019018290420.0009042791153430 -0.000050226371419-0.000100430028501-0.000150610971371-0.000200769200158-0.0004512196372670.0014066564016450.0019090336879470.0024114109742490.0049232974057590.0074351838372680.009947070268778-

27、0.000701102241287-0.000950417028058-0.001199164013417-0.002434382472993-0.003655408270323-0.0048622433804000.0124589567002880.0249589567002880.0374589567002880.0499589567002880.0624589567002880.074958956700288-0.006054889775763-0.011778983079828-0.017151998201403-0.022174175468426-0.026845753781789-

28、0.0311669705745320.0874589567002880.0999589567002880.1124589567002880.1249589567002880.1374589567002880.149958956700288-0.035138061683373-0.038759261501758-0.042030803032876-0.044952917848809-0.047525835962155-0.0497497859783920.1624589567002880.1749589567002880.1874589567002880.1999589567002880.2124589567002880.224958956700288-0.05162