欧拉格式 改进欧拉格式 后退欧拉格式
MATLAB改进欧拉法与四阶龙格-库塔求解一阶常微分方程
姓名：樊元君学号：02 日期：一、实验目的掌握MATLAB语言、C/C++语言编写计算程序的方法、掌握改进欧拉法与四阶龙格-库塔求解一阶常微分方程的初值问题。
掌握使用MATLAB程序求解常微分方程问题的方法。
：二、实验内容1、分别写出改进欧拉法与四阶龙格-库塔求解的算法，编写程序上机调试出结果，要求所编程序适用于任何一阶常微分方程的数值解问题，即能解决这一类问题，而不是某一个问题。
实验中以下列数据验证程序的正确性。
求，步长h=。
*2、实验注意事项的精确解为，通过调整步长，观察结果的精度的变化^)三、程序流程图：●改进欧拉格式流程图：~|●四阶龙格库塔流程图：]四、源程序：●改进后欧拉格式程序源代码：function [] = GJOL(h,x0,y0,X,Y)format longh=input('h=');…x0=input('x0=');y0=input('y0=');disp('输入的范围是：');X=input('X=');Y=input('Y=');n=round((Y-X)/h);\i=1;x1=0;yp=0;yc=0;for i=1:1:nx1=x0+h;yp=y0+h*(-x0*(y0)^2);%yp=y0+h*(y0-2*x0/y0);%·yc=y0+h*(-x1*(yp)^2);%yc=y0+h*(yp-2*x1/yp);%y1=(yp+yc)/2;x0=x1;y0=y1;y=2/(1+x0^2);%y=sqrt(1+2*x0);%fprintf('结果=%.3f,%.8f,%.8f\n',x1,y1,y);:endend●四阶龙格库塔程序源代码：function [] = LGKT(h,x0,y0,X,Y)。
format longh=input('h=');x0=input('x0=');y0=input('y0=');disp('输入的范围是：');"X=input('X=');Y=input('Y=');n=round((Y-X)/h);i=1;x1=0;k1=0;k2=0;k3=0;k4=0;for i=1:1:n~x1=x0+h;k1=-x0*y0^2;%k1=y0-2*x0/y0;%k2=(-(x0+h/2)*(y0+h/2*k1)^2);%k2=(y0+h/2*k1)-2*(x0+h/2)/(y0+h/2*k1);% k3=(-(x0+h/2)*(y0+h/2*k2)^2);%k3=(y0+h/2*k2)-2*(x0+h/2)/(y0+h/2*k2);% k4=(-(x1)*(y0+h*k3)^2);%k4=(y0+h*k3)-2*(x1)/(y0+h*k3);%…y1=y0+h/6*(k1+2*k2+2*k3+k4);%y1=y0+h/6*(k1+2*k2+2*k3+k4);%x0=x1;y0=y1;y=2/(1+x0^2);%y=sqrt(1+2*x0);%fprintf('结果=%.3f,%.7f,%.7f\n',x1,y1,y);end·end*五、运行结果：改进欧拉格式结果：;}四阶龙格库塔结果：步长分别为：和时，不同结果显示验证了步长减少，对于精度的提高起到很大作用，有效数字位数明显增加。
四阶龙格库塔和向前欧拉方法和隐式欧拉法以及改进欧拉法教材
实验九欧拉方法信息与计算科学金融崔振威201002034031一、实验目的：1、掌握欧拉算法设计及程序实现二、实验内容：1、p364-9.2.4、p386-9.5.6三、实验要求：主程序：欧拉方法（前项）：function [x,y]=DEEuler(f,a,b,y0,n);%f:一阶常微分方程的一般表达式的右端函数%a:自变量的取值下限%b:自变量的取值上限%y0：函数的初值%n:积分的步数if nargin<5,n=50;endh=(b-a)/n;x(1)=a;y(1)=y0;for i=1:nx(i+1)=x(i)+h;y(i+1)=y(i)+h*feval(f,x(i),y(i));endformat short欧拉方法（后项）：function [x,y]=BAEuler(f,a,b,y0,n);%f:一阶常微分方程的一般表达式的右端函数%a:自变量的取值下限%b:自变量的取值上限%y0：函数的初值%n:积分的步数if nargin<5,n=50;endh=(b-a)/n;x(1)=a;y(1)=y0;for i=1:nx(i+1)=x(i)+h;y1=y(i)+h*feval(f,x(i),y(i));y(i+1)=y(i)+h*feval(f,x(i+1),y1);endformat short梯形算法：function [I,step,h2] = CombineTraprl(f,a,b,eps)%f 被积函数%a,b 积分上下限%eps 精度%I 积分结果%step 积分的子区间数if(nargin ==3)eps=1.0e-4;endn=1;h=(b-a)/2;I1=0;I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h;while abs(I2-I1)>epsn=n+1;h=(b-a)/n;I1=I2;I2=0;for i=0:n-1x=a+h*i;x1=x+h;I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+subs(sym(f),findsym(sym(f)),x1));endendI=I2;step=n;h2=(b-a)/n;改进欧拉方法：function [x,y]=MoEuler(f,a,b,y0,n);%f:一阶常微分方程的一般表达式的右端函数%a:自变量的取值下限%b:自变量的取值上限%y0：函数的初值%n:积分的步数if nargin<5,n=50;endh=(b-a)/n;x(1)=a;y(1)=y0; for i=1:nx(i+1)=x(i)+h;y1=y(i)+h*feval(f,x(i),y(i)); y2=y(i)+h*feval(f,x(i+1),y1); y(i+1)=(y1+y2)/2; endformat short四阶龙格-库塔法：function y = DELGKT4_lungkuta(f, h,a,b,y0,varvec) %f:一阶常微分方程的一般表达式的右端函数 %h:积分步长%a ：自变量的取值下限 %b:自变量的取值上限%varvec ：常微分方程的变量组 format long; N = (b-a)/h;y = zeros(N+1,1); y(1) = y0; x = a:h:b;var = findsym(f); for i=2:N+1K1 = Funval(f,varvec,[x(i-1) y(i-1)]);%Funval 为程序所需要的函数 K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K2*h/2]); K4 = Funval(f,varvec,[x(i-1)+h y(i-1)+h*K3]); y(i) = y(i-1)+h*(K1+2*K2+2*K3+K4)/6; endformat short;p364-9.2.4欧拉方法（前项）：1、22)(,1)0(,22'+-+-==-=-t t e t y y y t y t解:执行20步时：编写函数文件doty.m 如下： function f=doty(x,y); f=x^2-y;在Matlab 命令窗口输入：>> [x1,y1]=DEEuler('doty',0,2,1,20)可得到结果：x1 =0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.70000.8000 0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.50001.6000 1.7000 1.8000 1.90002.0000y1 =1.0000 0.9000 0.8110 0.7339 0.6695 0.6186 0.5817 0.55950.5526 0.5613 0.5862 0.6276 0.6858 0.7612 0.8541 0.96471.0932 1.2399 1.4049 1.5884 1.7906在Matlab 命令窗口输入：>> y3=-exp(-x1)+x1.^2-2*x1+2求得解析解：y3 =1.0000 0.9052 0.8213 0.7492 0.6897 0.6435 0.6112 0.59340.5907 0.6034 0.6321 0.6771 0.7388 0.8175 0.9134 1.02691.1581 1.3073 1.4747 1.6604 1.8647输入：>> plot(x1,y1,'o');hold on>> plot(x1,y3,'*');hold on可得近似值与解析解的图像比较：执行40步时：在Matlab 命令窗口输入：>> [x1,y1]=DEEuler('doty',0,2,1,40)可得到结果：x1 =0 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 0.6000 0.6500 0.7000 0.75000.8000 0.8500 0.9000 0.9500 1.0000 1.0500 1.1000 1.15001.2000 1.2500 1.3000 1.3500 1.4000 1.4500 1.5000 1.55001.6000 1.6500 1.7000 1.7500 1.8000 1.8500 1.9000 1.95002.0000y1 =1.0000 0.9500 0.9026 0.8580 0.8162 0.7774 0.7417 0.7091 0.6798 0.6538 0.6312 0.6121 0.5967 0.5848 0.5767 0.5724 0.5719 0.5753 0.5826 0.5940 0.6094 0.6290 0.6526 0.68050.7126 0.7490 0.7897 0.8347 0.8841 0.9379 0.9961 1.05881.1260 1.1977 1.2739 1.3547 1.4401 1.5301 1.6247 1.7240 1.8279在Matlab 命令窗口输入：>> y3=-exp(-x1)+x1.^2-2*x1+2求得解析解：y3 =1.0000 0.9513 0.9052 0.8618 0.8213 0.7837 0.7492 0.7178 0.6897 0.6649 0.6435 0.6256 0.6112 0.6005 0.5934 0.5901 0.5907 0.5951 0.6034 0.6158 0.6321 0.6526 0.6771 0.70590.7388 0.7760 0.8175 0.8633 0.9134 0.9679 1.0269 1.09031.1581 1.2305 1.3073 1.3887 1.4747 1.5653 1.6604 1.7602 1.8647输入：>> plot(x1,y1,'o');hold on>> plot(x1,y3,'*');hold on可得近似值与解析解的图像比较：从上面结果可以看出，执行40步时近似解的值要接近于解析解，误差更小，结果更精确。
第五章：常微分方程数值解法第一节欧拉法
改进欧拉法公式
改进欧拉法公式
拓扑学中欧拉公式的证明
用数学归纳法证明
( 1)当 r= 2时，由表明 1,这两个区域可以想象为以赤道为边界的两个半球面，
赤道上存有两个“顶点” 将赤道分为两条“边界”，即为 r= 2，v= 2，e= 2；于是 r+
v- e= 2,欧拉定理设立.。
( 2)设r= m(m≥ 2)时欧拉定理成立，下面证明 r= m+ 1时欧拉定理也成立。
由表明 2，我们在 r= m+ 1的地图上自由选择一个区域 x ,则 x 必存有与它如此相连的区域 y ，使在换成 x 和 y 之间的唯一一条边界后，地图上只有 m 个区域了；在
换成 x 和 y 之间的边界后，若原该边界两端的顶点现在都还是 3条或 3条以上边界
的顶点，则该顶点留存，同时其他的边界数维持不变;若原该边界一端的或两端的顶
点现在沦为 2条边界的顶点，则换成该顶点 ,该顶点两边的两条边界便沦为一条边界。
于就是 ,在换成 x 和 y之间的唯一一条边界时只有三种情况：
①减少一个区域和一条边界；
②增加一个区域、一个顶点和两条边界；
③减少一个区域、两个顶点和三条边界；
即为在换成 x 和 y 之间的边界时 ,不论何种情况都必定存有“增加的区域数 + 增
加的顶点数 = 增加的边界数”我们将上述过程反过来 (即将 x 和 y之间换成的边界又
照曝光原样图画上 ) ，就又沦为 r= m+ 1的地图了，在这一过程中必然就是“减少的
区域数 + 减少的顶点数 = 减少的边界数”。
因此，若r= m (m≥2)时欧拉定理成立，则 r= m+ 1时欧拉定理也成立.。
由 ( 1)和 ( 2)所述 ,对于任何正整数r≥2,欧拉定理设立。
复变函数中欧拉公式的证明。
欧拉格式 改进欧拉格式 后退欧拉格式
functionf=doty(x,y)
f=y-2*x/y;
在matlab命令窗口输入：
[x,y]=euler('doty',0,1,1,10);
y1=sqrt(1+2.*x); % 对应的准确解
disp(' x y y1')
disp([x',y',y1'])
(2)、后退的Euler格式
e3=[0 0.0005 0.0009 0.0013 0.0017 0.0022 0.0027 0.0033 0.0040 0.0048 0.0058];
plot(x,abs(e1),'*');
hold on
plot(x,abs(e2),'b');
hold on
plot(x,abs(e3),'r-');
5、实验总结
通过编程实现了常微分方程初值问题数值解法中的欧拉方法及其后退、改进的算法，并比较了其数值解与精确解之间的误差。可以看出后退的欧拉方法得到的数值解精确度较差，而改进的欧拉方法得到的结果则相对较好。
6、教师评语及评分
hold on
4、实验结果与分析
（1）、Euler格式
计算常微分方程的结果为：
>> [x,y]=euler('doty',0,1,1,10);
y1=sqrt(1+2.*x); % 对应的准确解
disp(' x y y1')
disp([x',y',y1'])
x y y1
0 1.0000 1.0000
0.1000 1.1000 1.0954
Euler格式和RungeKutta格式.ppt
讨论常微分方程（ordinary differential equation）的定解问题，这类问题
的最简单型是如下一阶方程的初值问题：
y f ( x, y), a x b
y( x0 )
y0 ,
y .
(*)
由于数值解是找精确解 y( x) 的近似值，因此，总假设方程的解 y( x)
可用拟合方法求该组数据 ( xi , yi ) 的近似曲线
➢得到数值解有两个基本途径：
⑴ 把近似解表示成有限个独立函数之和，例如：截断的幂（Taylor）
级数或正交函数展开式中的前几项. 涉及到计算高阶导，尽管可用“数 值微分”技术，但得到的公式太长、太复杂！通常比较适用于手算.
⑵ 离散化方法（也称为差分方法），它提供了用当前节点上的或前几
存在且唯一，并具有充分的光滑性！
定理 设函数f(x,y)在区域 D : a x b, y R 上连续,且在区域D内关于y
满足Lipschitz (李普希兹)条件，即存在正数L，使得对所有的 x [a, b]
和任何 y1, y均2有
f ( x, y1 ) f ( x, y2 ) L y1 y2 ,
ba n
②从初始条件 y( x0 ) y0 出发，按照节点的排序，依次逐个计算 yi 的
值，称为步进法. 一般有两种类型：单步法、多步法.
➢为了考察数值解是否具有使用价值，必须解决的基本问题：
①当步长h取得充分小时，所得的近似解yn能否以足够的精度逼近初值
问题的精确解y(xn).这就是收敛问题。即当h 0时, yn y(xn) ?
再从 ( x1, y1) 出发，以 f ( x1, y1) 为斜率作直线 y y1 f ( x1, y1)( x x1)
偏微分方程数值解例题答案
yyy[y例11110.1[1(101)]0.9,10.1[0.9(10.10.9)]0.9019,1(0.90.9019)0.900952p c y y y ì=-´´+´=ïï=-´´+´=íïï=+=î20.900950.1[0.90095(10.10.90095)]0.80274,0.900950.1[0.80274(10.20.80274)]0.80779,1(0.802740.80779)0.805262p c y y yì=-´´+´=ïï=-´´+´=íïï=+=î 这样继续计算下去，其结果列于表9.1. 表9.1 Euler 方法方法改进的Euler 方法方法准确值准确值n xn yny)(n x y0.1 0.9000000 0.9009500 0.9006235 0.2 0.8019000 0.8052632 0.8046311 0.3 0.7088491 0.7153279 0.7144298 0.4 0.6228902 0.6325651 0.6314529 0.5 0.5450815 0.5576153 0.5563460 0.6 0.4757177 0.4905510 0.4891800 0.7 0.4145675 0.4310681 0.4296445 0.8 0.3610801 0.3786397 0.3772045 0.9 0.3145418 0.3326278 0.3312129 1.0 0.2741833 0.2923593 0.2909884 从表9.1可以看出,Euler 方法的计算结果只有2位有效数字,而改进的Euler 方法确有3位有效数字,这表明改进的Euler 方法的精度比Euler 方法高. 例2 试用Euler 方法、改进的Euler 方法及四阶经典R-K 方法在不同步长下计算初值问题ïîïíì=££+-=1)0(,10),1(d d y x xy y xy 在0.2、0.4、0.8、1.0处的近似值，并比较它们的数值结果. 解 对上述三种方法,每执行一步所需计算)1(),(xy y y x f +-=的次数分别为1、2、4。
用改进欧拉法和梯形法解初值问题
在数值分析领域，求解初值问题是一个非常重要的课题。
常见的数值方法有欧拉法和梯形法，但它们都存在一定的局限性。
本文将深入探讨如何改进欧拉法和梯形法，以更高效、精确地解决初值问题。
让我们回顾一下初值问题的基本概念。
初值问题是指，在已知微分方程初值条件的情况下，求解微分方程的数值解。
通常情况下，微分方程很难直接求解，因此需要借助数值方法来逼近微分方程的解。
欧拉法和梯形法是最常见的数值方法之一，它们都是通过离散化微分方程来逼近微分方程的解。
接下来，让我们讨论欧拉法的改进方法。
欧拉法是一种一阶精度的数值方法，它的局限性在于步长选择较大时，数值解会出现较大的误差。
为了改进欧拉法，可以考虑使用改进的欧拉法，如改进的Euler-Cauchy法，它通过使用更复杂的插值公式来提高精度。
也可以考虑使用四阶Runge-Kutta法等更高阶的方法来提高数值解的精度和稳定性。
另让我们分析梯形法的改进方法。
梯形法是一种二阶精度的数值方法，它对微分方程积分过程进行了更加精确的近似。
但梯形法在处理刚性微分方程时会出现数值不稳定的情况。
为了改进梯形法，可以考虑使用隐式的Runge-Kutta法，它可以更好地处理刚性微分方程，提高数值解的稳定性和精度。
总结来说，改进欧拉法和梯形法是为了解决数值解精度不高、数值不稳定等问题。
在实际应用中，根据不同的初值问题特性，选择合适的数值方法很关键。
只有综合考虑数值解精度、稳定性等因素，才能更好地求解初值问题。
从我个人的观点来看，改进欧拉法和梯形法是数值分析领域的一个重要研究方向。
随着科学技术的不断发展，解决初值问题需要更高效、精确的数值方法。
改进欧拉法和梯形法对于提高数值解的精度和稳定性具有重要意义。
在本文中，我们深入探讨了如何改进欧拉法和梯形法，以更高效、精确地解决初值问题。
通过对不同的改进方法进行分析，我们可以更好地理解数值方法的优缺点，为实际问题的求解提供参考。
希望本文能够对初值问题的数值解求解提供一定的帮助，也希望读者可以从中获得一定的启发和思考。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
>> x=0:0.1:1;
e1=[0 0.0046 0.0086 0.0125 0.0166 0.0209 0.0257 0.0311 0.0373 0.0445 0.0527];
e2=[0 -0.0046 -0.0089 -0.0132 -0.0179 -0.0232 -0.0293 -0.0364 -0.0449 -0.0551 -0.0674];
1.0000 1.7848 1.7321
(2)、后退的Euler格式
计算常微分方程的结果为：
>> [x,y]=backeuler('doty',0,1,1,10);
y1=sqrt(1+2.*x); % 对应的准确解
disp(' x y y1')
disp([x',y',y1'])
x y y1
0 1.0000 1.0000
0.8000 1.6165 1.6125
0.9000 1.6782 1.6733
1.0000 1.7379 1.7321
（4）、三者的误差图
注：‘*’表示Euler格式的数值解与准确解的误差，‘b’表示后退的Euler格式的数值解与准确解的误差，‘r-’表示改进的Euler格式的数值解与准确解的误差。
5、实验总结
通过编程实现了常微分方程初值问题数值解法中的欧拉方法及其后退、改进的算法，并比较了其数值解与精确解之间的误差。可以看出后退的欧拉方法得到的数值解精确度较差，而改进的欧拉方法得到的结果则相对较好。
6、教师评语及评分
（4）、根据结果，绘制Euler格式、后退的Euler格式和改进的Euler格式三者的误差图
根据表1可得：
Euler格式的数值解与准确解之间的误差e1：
e1=y-y1=[0 0.0046 0.0086 0.0125 0.0166 0.0209 0.0257 0.0311 0.0373 0.04450.0527];
1.4832
0.7
1.5803
1.5128
1.5525
1.5492
0.8
1.6498
1.5676
1.6165
1.6125
0.9
1.7178
1.6183
1.6782
1.6733
1.0
1.7848
1.6647
1.7379
1.7321
3、详细设计
(1)、Euler格式
欧拉格式的matlab实现程序（euler.m）如下：
微分方程数值解实验报告
实验序号：1日期：2013年10月15日
班级
10应数A班
姓名
黄达
学号
201005050148
实验名称
Euler格式
实验所用软件及版本
Maltab2008
1、实验目的
进一步理解Euler格式、后退的Euler格式、改进的Euler格式的设计思路和算法流程，培养动手实践能力和分析能力。
function[x,y]=neweuler(fun,x0,xfinal,y0,n)
ifnargin<5,n=10;
end
h=(xfinal-x0)/n;
x(1)=x0;y(1)=y0;
fori=1:n
x(i+1)=x(i)+h;
z0=y(i)+h*feval(fun,x(i),y(i));
y(i+1)=y(i)+h/2*(feval(fun,x(i),y(i))+feval(fun,x(i+1),z0));
end
编写一个doty.m文件，如下：
functionf=doty(x,y)
f=y-2*x/y;
在matlab命令窗口输入：
[x,y]=neweuler('doty',0,1,1,10);
y1=sqrt(1+2.*x); % 对应的准确解
disp(' x y y1')
disp([x',y',y1'])
f=y-2*x/y;
在matlab命令窗口输入：
[x,y]=backeuler('doty',0,1,1,10);
y1=sqrt(1+2.*x); % 对应的准确解
disp(' x y y1')
disp([x',y',y1'])
(3)、改进的Euler格式
改进的Euler格式的matlab实现程序（neweuler.m）如下：
后退的Euler格式的matlab实现程序（backeuler.m）如下：
function[x,y]=backeuler(fun,x0,xfinal,y0,n)
ifnargin<5,n=10;
end
h=(xfinal-x0)/n;
x(1)=x0;y(1)=y0;
fori=1:n
x(i+1)=x(i)+h;
1.1000
1.0909
1.0959
1.0954
0.2
1.1918
1.1743
1.1841
1.1832
0.3
1.2774
1.2517
1.2662
1.2649
0.4
1.3582
1.3237
1.3434
1.3416
0.5
1.4351
1.3910
1.4164
1.4142
0.6
1.5090
1.4540
1.4860
0.1000 1.0909 1.0954
0.2000 1.1743 1.1832
0.3000 1.2517 1.2649
0.4000 1.3237 1.3416
0.5000 1.3910 1.4142
0.6000 1.4540 1.4832
0.7000 1.5128 1.5492
0.8000 1.5676 1.6125
2、实验内容
编写Euler格式、后退的Euler格式、改进的Euler格式的程序代码，用于计算下列常微分方程，将计算结果列于表1，并绘制三者的误差图，给出相应的结论。
准确解 ，步长取 。
表1. Euler格式的数值结果（保留5位有效数字）
Euler格式
后退Euler
格式
改进的Euler格式
准确解
0.1
e3=[0 0.0005 0.0009 0.0013 0.0017 0.0022 0.0027 0.0033 0.0040 0.0048 0.0058];
plot(x,abs(e1),'*');
hold on
plot(x,abs(e2),'b');
hold on
plot(x,abs(e3),'r-');
0.2000 1.1918 1.1832
0.3000 1.2774 1.2649
0.4000 1.3582 1.3416
0.5000 1.4351 1.4142
0.6000 1.5090 1.4832
0.7000 1.5803 1.5492
0.8000 1.6498 1.6125
0.9000 1.7178 1.6733
z0=y(i)+h*feval(fun,x(i),y(i));
fork=1:3
z1=y(i)+h*feval(fun,x(i+1),z0);
ifabs(z1-z0)<1e-3
break;
end
z0=z1;
end
y(i+1)=z1;
end
编写一个doty.m文件，如下：
functionf=doty(x,y)
后退的Euler格式的数值解与准确解之间的误差e2：
e2=y-y1=[0 -0.0046 -0.0089 -0.0132 -0.0179 -0.0232 -0.0293 -0.0364 -0.0449 -0.0551-0.0674];
改进的Euler格式的数值解与准确解之间的误差e2：
e3=y-y1=[0 0.0005 0.0009 0.0013 0.0017 0.0022 0.0027 0.0033 0.0040 0.0048 0.0058].
0.9000 1.6183 1.6733
1.0000 1.6647 1.7321
（3）、改进的Euler格式
计算常微分方程的结果为：
>> [x,y]=neweuler('doty',0,1,1,10);
y1=sqrt(1+2.*x); % 对应的准确解
disp(' x y y1')
disp([x',y',y1'])
function[x,y]=euler(fun,x0,xfinal,y0,n)
ifnargin<5,n=10;
end
h=(xfinal-x0)/n;
x(1)=x0;y(1)=y0;
fori=1:n
x(i+1)=x(i)+h;
y(i+1)=y(i)+h*feval(fun,x(i),y(i));
end
x y y1
0 1.0000 1.0000
0.1000 1.0959 1.0954
0.2000 1.1841 1.1832
0.3000 1.2662 1.2649
0.4000 1.3434 1.3416
0.5000 1.4164 1.4142
0.6000 1.4860 1.4832