常微分方程的数值解法
专业班级：信息软件 姓名：吴中原 学号：120108010002 一、实验目的
1、熟悉各种初值问题的算法，编出算法程序；
2、明确各种算法的精度与所选步长有密切关系；通过计算更加了解各种 算法的优越性。
二、实验题目
1、根据初值问题数值算法，分别选择二个初值问题编程计算；
x 2、试分别取不同步长，考察某节点j处数值解的误差变化情况； 3、试用不同算法求解某初值问题，结果有何异常； 4、分析各个算法的优缺点。
三、实验原理与理论基础
（一） 欧拉法算法设计
对常微分方程初始问题
?dy?f(x,y) ?dx?? ?y(x0)?y0
(6-1)
(6-2)
用数值方法求解时，我们总是认为(6-1)、(6-2)的解存在且唯一。
欧拉法是解初值问题的最简单的数值方法。从(6-2)式由于y (x0) = y0已给定，因而可以算出y'(x0)?f(x0,y0)。
设x1 = h充分小，则近似地有：
y(x1)?y(x0)?y'(x0)?f(x0,y0)
h (6-3)
记 yi?y(xi) i?0,1,?,n 从而我们可以取
y1?y0?hf(x0,y0)
作为y(x1)的近似值。利用y1及f (x1, y1)又可以算出y(x2)的近似值：
y2?y1?hf(x1,y1)
一般地，在任意点xn?1??n?1?h处y(x)的近似值由下式给出
yn?1?yn?hf(xn,yn)
这就是欧拉法的计算公式，h称为步长。
(6-4)
（二）四阶龙格-库塔法算法设计：
??yi?1?yi?k1 欧拉公式可以改写为：?，它每一步计算f?x,y?的值一次，截
??k1?hf?xi,yi?断误差为o?h2?。
1?y?y??k1?k2?i?i?12?k1?hf?xi,yi?改进的欧拉公式可以改写为：? ，它每一步要计算?k?hf?x?h,y?k?ii1?2?f?x,y?的值两次，截断误差为o?h3?。
改进的欧拉方法之所以比欧拉方法具有更高的精度，是因为在每一步它都比欧拉方法多计算了一次f?x,y?的值。因此，要进一步提高精度，可以考虑在每一步增加计算f?x,y?的次数。
如果考虑在每一步计算f?x,y?的值四次，则可以推得如下公式：
?1y?y??k1?2k2?2k3?k4??i?1i6??k1?hf?xi,yi? ?h1??? k?hfx?,y?k1??2i?i22????h1???k3?hf?xi?,yi?k2?22????k?hf?x?h,y?k?ii3?4此公式称为标准四阶龙格-库塔(Runge-Kutta)公式，它的截断误差为o?h5?。虽然用龙格-库塔方法每一步需要四次调用f，计算量较改进的欧拉方法大一倍，这里由于龙格-库塔方法的步长增大了一倍，因而两种方法总的计算量相同，但龙格-库塔方法精确度更高。所以龙格-库塔公式兼顾了精度和计算工作量的较为理想的公式，在实际计算中最为常用。
四、实验内容
（一）问题重述：
科学计算中经常遇到微分方程（组）初值问题，需要利用Euler法，改进Euler 法，Rung-Kutta方法求其数值解，诸如以下问题：
??? y??4xy ? xy（1）?? y?0? ? 0 0 < x?1
分别取h=0.1,0.2,0.4时数值解。
2 初值问题的精确解y?4?5e?x。
（2）用r=3的Adams显式和预 - 校式求解
? y ??x2?y2?? ??y??1??0取步长h=0.1，用四阶标准 ?1?x?0 R-K
方法求值。 （3）用改进Euler法或四阶标准R-K方法求解
?? y1?? y2 y1??0???1? y???y1 y?22?0?? 0
? y?3??y3 y3?0??1
0?x?1
取步长h?0.01,计算y?0.05?, y?0.10?, y?0.15?数值解，y1?0.15???0.9880787, y2?0.15??0.1493359, y3?0.15??0.8613125 （4）利用四阶标准R- K方法求二阶方程初值问题的数值解
??y???3y??2y?0 （I）?y?0??0,y??0??1 0?x?1 h?0.02
??y???0.1(1?y2)y??y?0 (II) ?y?0??1, y??0??0 0?x?1, h?0.1
??y?y?? (III) ?ex?1?y?0??1,y??0??0 0?x?2, h?0.1
??y???siny?0 (IV) ?y?0??1,y??0??0 0?x?4, h?0.2
(二)实验代码： 1、欧拉法程序
function y=Euler(a,b,M,y0)
%a=1,b=2,M=10,f=t*y^(1/3),y0=1; h=(b-a)/M;
t=zeros(1,M+1); t=a:h:b;
y=zeros(1,M+1); yy=zeros(1,M+1); y(1)=y0; for k=1:M
y(k+1)=y(k)+h*t(k)*y(k)^(1/3); end
yb=y(M+1);
参考结果