Yn+1=Yn+hf（Xn+th,Y（Xn+th））
令K=f（Xn+th,Y（Xn+th））为平均斜率。
从而有：
K1=f（Xn,Yn）;
K2=f（Xn+h/2,Yn+（1/2）hK1）;
Yn+1=Yn+h*（c1K1+c2K2）;
上式即为二阶Runge-Kutta方法，其中对c1,c2进行一定的取值，可得到更加简单的结果。
4、四阶Runge-Kutta方法
如果我们将上述方法加以推广，考虑初值问题
（1）的如下离散变量法：
Yn+1=Yn+hm（Xn,Yn,h）,
其中
m（Xn,Yn,h）=∑cr*Kr,
K1=f（Xn,Yn）
Kr=f（Xn+arh,Yn+∑（j=1,…,n-1）brj*Kj）
生活中常用4级四阶Runge-Kutta方法：
Yn=1=Yn+（h/6）*（K1+2K2+2K3+K4）
其中
K1=f（Xn,Yn）
K2=f（Xn+（1/2）*h，Yn+（1/2）*h*K1）
K3=f（Xn+（1/2）*h，Yn+（1/2）*h*K2）
K4=f（Xn+h,Yn+h*K3）
X0=a,Xn=a+nh,h=（b-a）/N,n=1,2,…,N.
此方法应用最广泛，最早。
三、实验程序：
1、euler方法：
建立m文件
function[x,y]=euler（fun,x0,xfinal,y0,N）;
h=（xfinal-x0）/N;
x
（1）=x0;
y
（1）=y0;
fori=1:
N
y（i+1）=y（i）+h*feval（fun,x（i）,y（i））;
x（i+1）=x（i）+i*h;
end
2、改进euler方法
function[x,y]=eulerchange（fun,x0,xfinal,y0,N）;
h=（xfinal-x0）/N;
x
（1）=x0;y
（1）=y0;
fori=1:
N
x（i+1）=x（i）+i*h;
y1=y（i）+h*feval（fun,x（i）,y（i））;
y2=y（i）+h*feval（fun,x（i+1）,y1）
y（i+1）=（y1+y2）/2;
end
3、Runge-Kutta方法
function[x,y]=rk2（fun,x0,xfinal,y0,N）;
h=（xfinal-x0）/N;
a=1/2;
fori=1:
N
K1=feval（fun,x（i）,y（i））;
K2=feval（fun,（x（i）+a*h）,（y（i）+a*h*K1））;
c1=（2*a-1）/（2*a）;c2=1/（2*a）;
y（i+1）=y（i）+h*（c1*K1+c2*K2）;
end
4、四阶Runge-Kutta方法
function[x,y]=Euler1（a,b,n,w,）
h=（b-a）/n;
x=a;
y=w;
fori=1:
n
k1=-y（i）^2+x（i）^2+1;
k2=-（y（i）+（h*K1）/2）^2+（x（i）+h/2）^2+1;
K3=-（y（i）+（h*k2）/2）^2+（x（i）+h/2）^2+1;
k4=-（y（i）+h*k3）^2+（x（i）+h）^2+1;
y（i+1）=y（i）+（k1+2*k2+2*k3+k4）*h/6;
x（i+1）=x（i）+i*h;
end
四、实验题目：
对初值问题
（1）取
，分别用Euler方法、改进Euler方法、二阶Runge-Kutta方法、四阶Runge-Kutta方法求解，并与精确解进行比较。
（2）分析上述方法的计算量和计算精度，体会四阶Runge-Kutta方法的优点
解：
（1）建立m文件’doty’:
functionf=doty（x,y）
f=-（y^2）+x^2+1;
>>[x,y]=euler（'doty',0,1,1,20）
x=
Columns1through10
00.05000.15000.30000.50000.75001.05001.40001.80002.2500
Columns11through20
2.75003.30003.90004.55005.25006.00006.80007.65008.55009.5000
Column21
10.5000
y=
Columns1through10
1.00001.00001.00011.00121.00561.01761.04391.09451.18261.3247
Columns11through20
1.54011.84962.27312.82523.51134.32295.23856.22847.26498.3311
Column21
9.4232
[x,y]=eulerchange（'doty',0,1,1,20）
y2=1.0001y2=1.0012y2=1.0050y2=1.0151y2=1.0372y2=1.0794y2=1.1530y2=1.2720y2=1.4529y2=1.7133y2=2.0695y2=2.5339y2=3.1118y2=3.7997y2=4.5861y2=5.4553y2=6.3914y2=7.3822y2=8.4198y2=9.5004
x=
Columns1through10
00.05000.15000.30000.50000.75001.05001.40001.80002.2500
Columns11through20
2.75003.30003.90004.55005.25006.00006.80007.65008.55009.5000
Column21
10.5000
y=
Columns1through10
1.00001.00011.00071.00341.01131.02991.06721.13421.24491.4162
Columns11through20
1.66612.01252.46923.04333.73254.52515.40336.34857.34588.3863
Column21
9.4664
>>[x,y]=rk2（‘doty’,0,1,1,20）
[x,y]=rk2（'doty',0,1,1,20）
?
?
?
Undefinedfunctionorvariable"x".
Errorin==>rk2at5
K1=feval（fun,x（i）,y（i））;
（程序存在问题，无法运行）
>>[x,y]=Euler1（0,1,20,1）
?
?
?
Error:
File:
Euler1.mLine:
2Column:
31
Unbalancedorunexpectedparenthesisorbracket.
（2）Euler方法是一阶方法，计算比较简单，精确度不高；改进的Euler方法为二阶方法，精确度较原方法高；二阶
Runge-Kutta方法计算量较大，精确度更高，四阶的最
为有效精确，应用更广泛，但计算更复杂。
五、实验分析及心得：
结果分析
Eluer方法的总误差是O（h）,称其为一阶方法。
改进的Euler方法的总体截断误差是O（h），这是二r阶方法。
一个方法的总体截断误差若为O（h），则，称其为r阶方法。
一般来说，方法的总体截断误差阶数越高，其精度越高。
龙格库塔法是通过计算不同点上的函数值，并对这些函数值作线性组合，构造近似公式，再把近似公式与解的泰勒展开式进行比较，使前面的若干项相同，从而使近似公式达到一定阶数。
二阶Runge-Kutta法是用k1和k2的加权平均值来近似k经推导得到二阶Runge-Kutta法最近本的形式
四阶Runge-Kutta法是用Runge-Kutta法有k1和k2k3和k4的加权平均值来近似k
在计算精度上，四阶经典龙格-库塔方法的误差最小，改进欧拉方法其次，欧拉方法误差则比较大，所以四阶经典龙格-库塔方法得到最佳的精度。
而在计算量上面，相应地，很明显的四阶经典龙格-库塔方法也是最大，改进欧拉方法其次，欧拉方法计算量最小。
这样的结果，说明了运用以上三种方法时，其计算量的多少与精度的大小成正比。
我们在实际运用与操作中，可以根据实际情况，选择这3种方法中的其中一种最适合的，追求精度的话，可以使用四阶经典龙格-库塔方法；而改进的欧拉方法，在精度上和计算量上都表现得很出色，能够满足一般情况；而欧拉方法更主要的是适用于对
的估计上，相应的，精度则有所欠缺。
个人心得体会
在整个课程设计完成的过程中，从始至终按照理论结合实践的思想，按照计划稳步进行。
课程设计内容清晰，理论详细，编写程序简单，达到了课程设计的目的及要求。
通过此次课程设计的过程，自身也收获了很多，不仅更深入的学习了理论知识，而且能够将书本上的理论知识结合程序运算实际应用，提高了动手实践能力。
同时也学习了很多与课程设计有关的计算机操作，尤其能够熟练应用Matlab软件。
通过“数值计算方法”和Matlab的学习，能够有效的将二者相结合解决一些数学上的实际应用问题。
这次课程设计获益良多相信不仅对于计算方法这门学科，对以后的工作学习都有极大的帮助。
实际成绩评定
签章
年月日
指导教师意见：
升级会员 