下载文档原格式
改进的欧拉公式与精确解的变化规律改进的欧拉公式是最常用的数值解法之一，它通过近似求解微分方程来得到数值解。
与精确解相比，改进的欧拉公式是通过将微分方程的导数从一个点近似为两个点的斜率来计算下一个点的数值解。
改进的欧拉公式的变化规律是随着步长的减小，数值解会更接近精确解。
这是因为当步长越小时，近似的斜率越接近真实的导数值，从而得到的数值解也更准确。
具体来说，改进的欧拉公式的变化规律可以描述为以下几点：1. 当步长减小时，数值解的误差也减小。
这意味着数值解更接近精确解。
2. 当步长趋近于零时，数值解逼近精确解。
这是因为在这种情况下，近似的斜率越来越接近真实的导数值，从而得到的数值解趋近于精确解。
3. 当步长增大时，数值解的误差也增大。
这是因为在这种情况下，近似的斜率与真实的导数值之间的差异会增大，导致数值解与精确解之间的差异也增大。
总之，改进的欧拉公式是一种数值解法，它可以通过近似求解微分方程来得到数值解。
随着步长的减小，数值解会更接近精确解。
在步长趋近于零的情况下，数值解逼近精确解。
当步长增大时，数值解的误差也增大。
进一步说明，改进的欧拉公式是欧拉公式的改进版，通过将微分方程的导数从一个点近似为两个点的斜率来提高数值解的准确性。
改进的欧拉公式可以写为以下形式：y_{n+1} = y_n + h \cdot \frac{f(x_n, y_n) + f(x_{n+1},y_{n+1})}{2}其中，y_n 是精确解在离散点 x_n 处的近似值，h 是步长，f(x_n, y_n) 是微分方程的导数。
改进的欧拉公式的准确度比欧拉公式更高，是因为它通过使用两个点的斜率的平均值来更准确地近似导数值。
改进的欧拉公式的变化规律可以归结为以下几点：1. 当步长 h 减小时，数值解的准确性提高。
这是因为较小的步长使得近似的斜率更接近真实的导数值，从而得到更精确的数值解。
2. 当步长 h 增大时，数值解的准确性降低。
这是因为较大的步长导致近似的斜率与真实的导数值之间的误差增加，从而导致数值解的误差变大。
实验九欧拉方法信息与计算科学金融崔振威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步时近似解的值要接近于解析解，误差更小，结果更精确。
分别利用欧拉法和改进欧拉法求解微分方程组的数值解欧拉法（Euler’s Method）和改进欧拉法（Improved Euler’s Method），是求解常微分方程数值解的两种常用方法。
它们都属于一阶精度的显式迭代算法。
首先，我们来介绍一下欧拉法。
欧拉法是一种简单的数值求解算法，它基于微分方程的定义，将微分方程转化为差分方程。
考虑一个一阶常微分方程 dy/dx = f(x, y)，并给定初始条件 y(x0)= y0，我们希望求解在给定区间 [x0, xn] 上方程的数值解。
首先，我们将区间 [x0, xn] 平均分成 N 个小区间，每个小区间的长度为 h = (xn - x0) / N。
然后，我们可以使用以下的欧拉迭代公式计算数值解：y[i+1] = y[i] + h * f(x[i], y[i])其中，x[i] = x0 + i * h，y[i] 是在点 x[i] 处的数值解。
通过不断迭代上述公式，我们可以获得[x0, xn] 上微分方程的数值解。
欧拉法的优点在于简单易懂，计算速度较快。
然而，欧拉法的缺点是精度较低，误差随着步长h 的增大而增大。
为了提高精度，我们可以使用改进欧拉法。
改进欧拉法，也称为龙格–库塔算法（Runge-Kutta Method）或四阶龙格–库塔方法，是一种基于欧拉法的改进算法。
改进欧拉法使用了更多的近似取值，以减小误差。
与欧拉法类似，我们将区间 [x0, xn] 平均分成 N 个小区间，每个小区间的长度为 h = (xn - x0) / N。
然后，我们可以使用以下的公式计算数值解：k1 = h * f(x[i], y[i])k2 = h * f(x[i] + h/2, y[i] + k1/2)y[i+1] = y[i] + k2其中，k1 和 k2 是计算过程中的辅助变量。
通过不断迭代上述公式，我们可以获得 [x0, xn] 上微分方程的数值解。
改进欧拉法相对于欧拉法而言，计算精度更高。
向前欧拉法,向后欧拉法与改进欧拉法求解微分方程
向前欧拉法、向后欧拉法和改进欧拉法是求解微分方程的常用数值方法。
这些方法都是基于欧拉公式，即将微分方程中的导数用差分代替，从而将微分方程转化为差分方程，进而用数值方法求解。
向前欧拉法是一种简单的数值方法，它利用当前时刻的导数来估计下一时刻的解。
具体来说，假设微分方程为dy/dt=f(y,t)，则向前欧拉法的迭代公式为：y_n+1=y_n+hf(y_n,t_n)，其中h为时间步长。
这个公式可以看作是在当前时刻上做一个切线，然后用这个切线的斜率来估计下一时刻的解。
向后欧拉法是一种更加精确的数值方法，它利用下一时刻的导数来估计当前时刻的解。
具体来说，向后欧拉法的迭代公式为：
y_n+1=y_n+hf(y_n+1,t_n+1)，其中h为时间步长。
这个公式可以看作是在下一时刻上做一个切线，然后用这个切线的斜率来估计当前时刻的解。
由于向后欧拉法需要解一个非线性方程，因此比向前欧拉法更加复杂。
改进欧拉法是向前欧拉法和向后欧拉法的结合，它利用当前时刻和下一时刻的导数来估计当前时刻的解。
具体来说，改进欧拉法的迭代公式为：y_n+1=y_n+(h/2)(f(y_n,t_n)+f(y_n+1,t_n+1))，其中h 为时间步长。
这个公式可以看作是在当前时刻和下一时刻上各做一个切线，然后将这两个切线的斜率取平均值来估计当前时刻的解。
改进欧拉法相对于向前欧拉法和向后欧拉法更加精确。
总的来说，向前欧拉法、向后欧拉法和改进欧拉法都是求解微分
方程的有力工具，使用时需要根据具体问题选择合适的方法。
用Euler法和改进的Euler法求u' -5u(0 < t w,，(0)=1的数值解，步长h=0.1, 0.05，并比较两个算法的精度。
解：1) 当步长h=0.1 时编写程序如下所示clf clear clc%直接求解微分方程y=dsolve( 'Dy=-5*y' , 'y(0)=1' , 't' )%Euler 法h=0.1;t=0:h:1; n=length(t);u=zeros(1,n); u(1)=1;zbu(1,1)=t(1); zbu(2,1)=u(1);for i=2:nf=-5*u(i-1); u(i)=u(i-1)+h*f;zbu(1,i)=t(i); zbu(2,i)=u(i);endzbu%改进的Euler 法v=zeros(1,n);v0=zeros(1,n); v(1)=1;zbv(1,1)=t(1); zbv(2,1)=v(1);for i=2:n f=-5*v(i-1); v0(i)=v(i-1)+h*f; v(i)=v(i-1)+h/2*(f-5*v0(i)); zbv(1,i)=t(i); zbv(2,i)=v(i);end zbvplot(t,u, 'r*' , 'markersize' ,10)hold on,plot(t,v, 'r.' , 'markersize' ,20)hold on,ezplot(y,[0,1])hold on,title( 'Euler 法和改进的Euler 法比较(h=0.1 )),grid onlegend( 'Euler 法’,'?改进的Euler 法','解析解')%解真值h=0.1;t=O:h:1;n=len gth(t);for i=1: ny(i)=1/exp(5*t(i)); %通过第一部分程序直接解得的解析解zby(1,i)=t(i);zby(2,i)=y(i);endzby我们可以得到计算后的结果图像如图一所示EulBr法和改进的Eutgr法比较)0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.9 09 1图1 Euler法和改进的Euler法比较（h=0.1）同时，我们得到Euler法，改进的Euler法和解析解的在各点处数值分别如下所示：t坐标0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 欧拉 1.0000 0.5000 0.2500 0.1250 0.0625 0.0313 0.0156 0.0078 0.0039 0.0020 0.0010 改进欧拉 1.0000 0.6250 0.3906 0.2441 0.1526 0.0954 0.0596 0.0373 0.0233 0.0146 0.0091 真值 1.0000 0.6065 0.3679 0.2231 0.1353 0.0821 0.0498 0.0302 0.0183 0.0111 0.0067 表1 Euler法和改进的Euler法在各点数值比较（h=0.1）为了比较Euler法和改进的Euler法的算法精度，在这里我们利用相对误差的概念进行评判。
P7731.利用改进欧拉方法计算下列初值问题，并画出近似解的草图：代码：%改进欧拉法function Euler(t0,y0,inv,h)n=round(inv(2)-inv(1))/h;t(1)=t0;y(1)=y0;for i=1:ny1(i+1)=y(i)+h*fun(t(i),y(i));t(i+1)=t(i)+h;y(i+1)=y(i)+1/2*h*(fun(t(i),y(i))+ fun(t(i+1),y1(i+1)))endplot(t,y,'*r')function y=fun(t,y);y=y+1;调用：Euler(0,3,[0,2],0.5)得到解析解：hold on;y=dsolve('Dy=y+1','(y(0)=3)','t');ezplot(y,[0,2])图像：代码：function Euler1(t0,y0,inv,h)n=round(inv(2)-inv(1))/h;t(1)=t0;y(1)=y0;for i=1:ny1(i+1)=y(i)+h*fun(t(i),y(i));t(i+1)=t(i)+h;y(i+1)=y(i)+1/2*h*(fun(t(i),y(i))+ fun(t(i+1),y1(i+1))) endplot(t,y,'*r')function y=fun(t,y);y=y^2-4*t;调用：Euler1(0,0.5,[0,2],0.2)图像：代码：function Euler2(t0,y0,inv,h)n=round(inv(2)-inv(1))/h;t(1)=t0;y(1)=y0;for i=1:ny1(i+1)=y(i)+h*fun(t(i),y(i));t(i+1)=t(i)+h;y(i+1)=y(i)+1/2*h*(fun(t(i),y(i))+ fun(t(i+1),y1(i+1))) endplot(t,y,'*r')function y=fun(t,y);y=(3-y)*(y+1);调用：Euler2(0,4,[0,5],1)得到解析解：hold on;y=dsolve('Dy=(3-y)*(y+1)','y(0)=4','t');ezplot(y)图像：代码：function Euler2(t0,y0,inv,h)n=round(inv(2)-inv(1))/h;t(1)=t0;y(1)=y0;for i=1:ny1(i+1)=y(i)+h*fun(t(i),y(i));t(i+1)=t(i)+h;y(i+1)=y(i)+1/2*h*(fun(t(i),y(i))+ fun(t(i+1),y1(i+1))) endplot(t,y,'*r')function y=fun(t,y);y=(3-y)*(y+1);调用：Euler2(0,4,[0,5],0.5)得到解析解：hold on;y=dsolve('Dy=(3-y)*(y+1)','y(0)=4','t');ezplot(y)图像：14.考虑满足初始条件（x(0)，y(0)）=(1,1)的下列方程组：选定时间步长 t=0.25，n=5.用改进欧拉方法求两个方程组的近似解；（1）代码：function Euler4(t0,int,n,h)t=t0;x(1)=int(1);y(1)=int(2);for i=1:nx1(i+1)=x(i)+h*xfun(t(i),x(i),y(i));y1(i+1)=y(i)+h*yfun(t(i),x(i),y(i));t(i+1)=t(i)+h;x(i+1)=x(i)+1/2*h*(xfun(t(i),x(i),y(i))+xfun(t(i+1),x1(i+1) ,y1(i+1)));y(i+1)=y(i)+1/2*h*(yfun(t(i),x(i),y(i))+yfun(t(i+1),x1(i+1) ,y1(i+1)));endplot(t,x,'o-r')hold onplot(t,y,'*-g')hold onplot(x,y)function x=xfun(t,x,y);x=y;function y=yfun(t,x,y);y=-2*x-3*y;调用函数：Euler4(0,[1,1],5,0.25)图像：（2）代码：function Euler5(t0,int,n,h)t=t0;x(1)=int(1);y(1)=int(2);for i=1:nx1(i+1)=x(i)+h*xfun(t(i),x(i),y(i));y1(i+1)=y(i)+h*yfun(t(i),x(i),y(i));t(i+1)=t(i)+h;x(i+1)=x(i)+1/2*h*(xfun(t(i),x(i),y(i))+xfun(t(i+1),x1(i+1) ,y1(i+1)));y(i+1)=y(i)+1/2*h*(yfun(t(i),x(i),y(i))+yfun(t(i+1),x1(i+1) ,y1(i+1)));endplot(t,x,'o-r')hold onplot(t,y,'*-g')hold onplot(x,y)function x=xfun(t,x,y);x=y+y^2;function y=yfun(t,x,y);y=-x+0.2*y-x*y+1.2*y^2;调用函数：Euler5(0,[1,1],5,0.25)图像：。
题目再现：
2. 当病人采取服用口服药或肌肉注射来治疗疾病时，药物虽然瞬间进入了体内，但它一般都集中与身体的某一部位，靠其表面与肌体接触而逐步被吸收。
假定身体系统是一个单房室系统，设t 时刻体内药物的总量为x(t)，则x(t)满足：
问题分析：运用欧拉公式求微分方程的数值解。
微分方程为：
11dx ,(0)0dt
k t k De kx x -=-= 1，欧拉折线法：
设h=0.1,即n=201时，1k 0.6=，k 0.2=，D=200.MATLAB 程序如下所示：
clear
f=sym('0.6*200*exp(-1*0.6*t)-0.2*x ');
a=0; b=20;
h=0.1;
n=(b-a)/h+1;
t=0; x=0;
szj=[t,x];
for i=1:n-1
11dx ,(0)0
dt k t k De kx x -=-=10011100 ,(,)()k t k k k k k k k
t x x x h f t x x h k De kx t t h -++==⎧⎪=+=+-⎨⎪=+⎩
x=x+h*subs(f,{'t','x'},{t,x});
t=t+h;
szj=[szj;t,x];
end
szj ;
x=dsolve(‘Dx=120*exp(-0.6*t)-0.2*x’,'x(0)=0','t')
T=[0:0.1:20];
X=subs(x,T);
plot(szj(:,1),szj(:,2),'or-',T,X,’b-’)
输出结果为：
解析解：x=
其中红线代表的是数值解，蓝线代表的是解析解。
可见，数值解的误差还是比较小的。
二、改进的欧拉法：有公式：
[]11100(,)(,)2
()i i i i i i h x x f t x f t x x x t +++=+
+= 与欧拉公式结合使用，有
(0)1(1)
()111(,)
[(,)(,)],0,1,2 (2)
i i i i k k i i i i i i x x hf t x h x x f t x f t x k +++++=+=++= 设h=0.1,即n=201时，1k 0.6=，k 0.2=，D=200.MATLAB 程序如下所示：
clear
f=sym('0.6*200*exp(-1*0.6*t)-0.2*x');
a=0; b=20;
h=0.1;
n=(b-a)/h+1;
t=0; x=0;
szj=[t,x];
for i=1:n-1
m=0;p=0;q=0;
m=x;
q=subs(f,{'t','x'},{t,m});
x=m+h*q;
if abs(p-x)>0.1
p=x;
x= m+h/2*(q+subs(f,{'t','x'},{t+h,x}));
end
t=t+h;
szj=[szj;t,x];
end
szj ;
x=dsolve('Dx=120*exp(-0.6*t)-0.2*x', 'x(0)=0','t')
T=[0:0.1:20];
X=subs(x,T);
plot(szj(:,1),szj(:,2),'or-',T,X,'b-')
可见，改进后的欧拉算法误差更小，极其精确，基本与原曲线重合。
