for(k=1:3)
z1=y(n)+h*feval(fun,x(n+1),z0);
if(abs(z1-z0)<1e-3)
break;
end
z0=z1;
end
y(n+1)=z1;%%Biblioteka 据向后欧拉公式计算y值end
T=[x',y']
改进欧拉代码如下：function[x,y]=Gaijineuler(f,x0,y0,xZ,h)
对于初值问题 稳定性的研究，易知其准确解为 ，假定初值经过扰动后变为 ，对于扰动后的解为 因此带来的扰动误差为 ，因此考虑 时 的值，它取决于 。易知，若 ，则原问题是稳定的；若 ，则原问题是不稳定的；若 ，则原问题是渐进稳定的。
实际遇到的大多数常微分方程初值问题都是稳定的，因此在后面的讨论数值解法时这常常是默认条件。
T=a:h:b;%计算各个分点
%f为一阶微分方程的函数；x0，y0为初始条件；xZ为取值范围的一个端点，h为区间步长
n=fix((xZ-x0)/h);%计算分点数
y(1)=y0;
x(1)=x0;
fori=1:n
x(i+1)=x0+i*h;
yp=y(i)+h*feval(f,x(i),y(i));
yc=y(i)+h*feval(f,x(i+1),yp);
% fun为一阶微分方程，x0，y0为初始条件，xN为取值范围的一个端点，N为区间个数x=zeros(1,N+1);
x=zeros(1,N+1);
x(1)=x0;
y(1)=y0;
h=(xN-x0)/N;%h为区间步长
for(n=1:N)
x(n+1)=x(n)+h;
z0=y(n)+h*feval(fun,x(n),y(n));
考虑一阶常微分方程的初值问题 ，如果存在实数 使得 则称f关于y满足利普希茨条件，L称为利普希茨常数。
对于常微分方程初值问题 ，考虑初值 的扰动是问题的解 发生偏差的情形。若 时 的偏差被控制在有界范围内，则称该初值问题是稳定的，否则该初值问题不稳定的。
特别地，若 时 的偏差收敛于零，则称该初值问题是渐进稳定的。
（3）比较欧拉法、改进欧拉法及龙格-库塔法，能够选择合适的方法进行问题的研究计算；
2实验内容：求微分方程 （欧拉法求解）
求微分方程 （改进欧拉法求解）
求微分方程 （龙格-库塔求解）
根据实验的结果进行分析，了解一般方法的的优缺点，稳定性，收敛性以及截断误差的分析，针对相应问题拿出有效方法得出最优的结果。
（4）对右端 用右矩形公式得 ，也叫隐式欧拉法。
误差分析：1.称 为计算 时的局部截断误差；
2.如果数值方法的局部截断误差为 ，那么称这种数值方法的阶数是p，其实p为非负整数。通常情况下，步长h越小，p越高，则局部截断误差越小，计算精
初泰勒展开有
则有 可见欧拉方法是一阶方法，精度不是很高。
2．改进欧拉方法：
x=zeros(1,N+1);
x=zeros(1,N+1);
x(1)=x0;
y(1)=y0;
h=(xN-x0)/N;%h为区间步长
for(n=1:N)
x(n+1)=x(n)+h;
y(n+1)=y(n)+h*feval(fun,x(n),y(n));%根据向前欧拉公式计算y值
end
T=[x',y']
（2）向后欧拉法：functiony=Euler2(fun,x0,y0,xN,N)
1．欧拉法：
依据：积分曲线上一点 的切线斜率等于函数值。
方法：推进法，初始点 出发，依照方向场在改点的方向推进到
向前欧拉法的得到:
(1)将 在 处泰勒展开 取h的线性部分，得
（2）将初值问题中得导数用向前差商来代替有 ，因此
（3）将 两边同时对x的区间 上积分 对右端 用左矩形公式得 ，此方法称向前欧拉法，也叫显示欧拉法。
二．相关背景知识介绍以及初值问题稳定性的研究：
在科学与工程问题中，常微分方程表述物理量的变化规律，应用非常广泛，比如，天体运动的轨迹，机器人控制，化学反应过程的描述和控制以及电路瞬态过程分析等。这些问题中要求解随时间变化的物理量，即位置函数 表示时间，而微分方程描述了未知函数与它的一阶或高阶导数之间的关系。
梯形公式：对右端 用梯形公式得+显然梯形公式是隐式公式。
改进欧拉公式：先用欧拉公式求的一个初步的近似值 ，成为预测值，预测值 的精度可能达不到要求，在用梯形公式将他校正一次，记为 ，这个结果成为校正值。
预测：
校正：
误差分析：记 为改进欧拉公式在 处的截断误差，
记 因此有
， 表示在 出的局部截断误差。由此得，梯形公式的局部截断误差为 ，因此改进欧拉的截断误差为 ，可见改进欧拉的方法是二阶方法，改进欧拉方法优于欧拉方法。
（一阶龙格-库塔）当r=1时，这就是欧拉法。
（二阶龙格-库塔）当r=2时， ，这就是改进欧拉法。
三阶和四阶龙格-库塔也只是在一般情况下得结果。
三阶
四阶
其局部截断误差是：
三．程序代码
欧拉法代码如下：
（1）向前欧拉法：
functiony=Euler1(fun,x0,y0,xN,N)
%fun为一阶微分方程，x0，y0为初始条件，xN为取值范围的一个端点，N为区间个数
y(i+1)=(yp+yc)/2;%根据改进欧拉公式计算结果
end
x=x';
y=y';
1．
3.龙格-库塔代码如下：
（三阶龙格-库塔）
functionR=Longgekuta3(f,a,b,aZ,h)
%a，b为端点，h为步长，aZ为初值
n=(b-a)/h;
T=zeros(1,n+1);%定义向量
Y=zeros(1,n+1);
常微分方程初值问题数值解法的比较
数值计算实践—课程设计报告
课题名称
常微分方程初值问题数值解法的比较
完成时间
2013-1-17
姓名
班级
学号
成绩
一. 实验目的及内容
1实验目的：（1）了解常微分方程初值问题的理论背景以及初值问题稳定性、收敛性的研究；
（2）熟练掌握欧拉法、改进欧拉法、龙格-库塔法以及截断误差分析；
3．龙格—库塔法:
根据拉格朗日微分中值定理， ，记 得到 ，这样，只要给出一种计算 的算法，就可以得到相应的计算公式。欧拉公式可以写为
改进欧拉公式可以写成 因此推出一般的推出广式
称为p阶龙格-库塔方法，简称p阶R-K方法。
因为 这里的 均为常数。
因为给定的系数不唯一，因此这里的常数有无穷多个解，下面是特殊情况下和一般情况下得结果