【微分方程数值解】改进Euler法的具体实现以及三种Euler法比较

实验目的：

用预估矫正实现改进Euler法的具体实现，并比较三种Euler法解决初值问题。

实验原理：

 实验步骤：

·改进Euler法的具体实现：

Step1：进行改进Euler的编写，我们由公式可知改进Euler法有预估和矫正的算法实现，首先进行对T、h、t、N参数的设置。

Step2：由实验原理中改进Euler法的公式

设置代码中的精确数字位数1e-6；

Step3：运行代码得出图像

·三种Euler法的比较：

Step1：将所写的改进Euler写入所给的example_Euler.m文件中。

Step2：将三种方法用不同的线条绘制于同一个图中与精确解进行比较。

Step3：运行代码得出图像

·分析与讨论

对于初值问题数值解法一般采取步进式求解，即通过前面求解的节点结果推出后续节点的值

若某算法的局部截断误差为则称该算法具有p阶精度。

·思考与回答

代码

改进Euler法的具体实现：

function gaijin
T = 1;
h = 0.01;
t = 0:h:T;
N = length(t)-1;
solu = exp(-5.0*t);
u0 = 1;
f = @f1;
u_modified = update_euler(f,u0,t,h,N);
figure (1)
plot(t,u_modified,'d',t,solu,'r')
legend('modeuler','Exact Solution')
end
function u = update_euler(f,u0,t,h,N)
u    = zeros(N+1,1);
u(1) = u0;
eps_in = 1e-6;
  K_in = 6;
for n = 1:N
    s1 = u(n);
    du = 1;
    k  = 1;
    while abs(du)>eps_in && k<K_in
        s2 = u(n)+h/2*(f(t(n+1),s1)+f(t(n),u(n)));
        du = s2-s1;
        s1 = s2;
         k = k+1;end
    u(n+1) = s2;
end
end
function f =f1(t,u)
f = -5*u;
End

三种Euler法的比较：

function example_Euler
T = 1;
h = 0.01;
t = 0:h:T;
N = length(t)-1;
solu = exp(-5.0*t);
u0 = 1;
f = @f1;
u_euler = euler(f,u0,t,h,N);
u_implicit = implicit_euler(f,u0,t,h,N);
u_modified = update_euler(f,u0,t,h,N);
figure (1)
plot(t,u_euler,'*r',t,u_implicit,'o',t,u_modified,'d',t, solu,'r')
legend('Euler','imeuler','modeuler','Exact Solution')
end
function u = euler(f,u0,t,h,N)
u    = zeros(N+1,1);
u(1) = u0;
for n = 1:N
    fn     = f(t(n),u(n));
    u(n+1) = u(n)+h*fn;
end
end
function u = implicit_euler(f,u0,t,h,N)
u    = zeros(N+1,1);
u(1) = u0;
eps_in = 1e-6;
  K_in = 6;
for n = 1:N
    s1 = u(n);
    du = 1;
    k  = 1;
    while abs(du)>eps_in && k<K_in
        s2 = u(n)+h*f(t(n+1),s1);
        du = s2-s1;
        s1 = s2;
         k = k+1;
    end
    u(n+1) = s2;
end
end
function u = update_euler(f,u0,t,h,N)
u    = zeros(N+1,1);
u(1) = u0;
eps_in = 1e-6;
  K_in = 6;
for n = 1:N
    s1 = u(n);
    du = 1;
    k  = 1;
    while abs(du)>eps_in && k<K_in
        s2 = u(n)+h/2*(f(t(n+1),s1)+f(t(n),u(n)));
        du = s2-s1;
        s1 = s2;
         k = k+1;end
    u(n+1) = s2;
end
end
function f =f1(t,u)
f = -5*u;
end