Euler法matlab表示

function example1_euler_method


T= 1;
h = 0.1;
% h = 0.05;
% h = 0.01;
t = 0:h:T;
N = length(t) - 1;
solu = exp(-5*t);

u0 = 1 ;
f = @f1;
u_euler = euler(f,u0,t,h,N);
u_implicit = implicit_euler(f,u0,t,h,N);
u_modified = modified_euler(f,u0,t,h,N);

figure(1)
plot(t,u_euler,'*b', t,u_implicit,'gs',t,u_modified,'s',t,solu,'r')
legend('Euler','Implicit-Euler','Modified-Euler','Exact-soln')
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 = modified_euler(f,u0,t,h,N)%改进的euler法
u = zeros(N+1,1);
u(1) = u0;
eps_in = 1e-6;
K_in = 6;
for n = 1:N
fn = f(t(n),u(n));
s1 = u(n);
du = 1;
k = 1;
while abs(du)>eps_in & k<K_in
    s2 = u(n) + 0.5*h*(fn + f(t(n+1),s1));
    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法中$y_n+1$可由$y_n$直接表示出，而隐Euler法和改进的Euler法中，$y_n+1$不可可由$y_n$直接表示出，使用不动点迭代法求解$y_n+1$.
代码中当精度小于1e-6或迭代满六次得出$y_n+1$。

因为h足够小，且$f^{{}^{\prime }}=-5$满足不动点定理条件。

以下为t与隐Euler法的图。

clc,clear
T= 1;
h = 0.1;
% h = 0.05;
% h = 0.01;
t = 0:h:T;
N = length(t) - 1;
solu = exp(-5*t);
u0 = 1 ;
f = @f1;

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


figure(1)
u_implicit = u;
plot(t,u_implicit,'gs',t,solu,'r')


function f = f1(t,u)
f =-5*u;
end

隐Euler法和改进的Euler法中，使用不动点迭代法求解$y_n+1$.