RungeKutta法 matlab实现

function chap1_RungeKutta_method
% test the Runge-Kutta method for ODE IVP
% foundate = '2015-3-8';
% chgedate = '2020-04-30';
% by Zhang, Xuping

u0 = 1;
T = 1;
h = 0.1;
N = T/h;
t = 0:h:T;
solu = exact1(t);

f = @f1;
u_euler = euler(f,u0,t,h,N);
u_RK = runge_kutta_44(f,u0,t,h,N);

figure(1)
plot(t,u_euler,'*b',t,u_RK,'o',t,solu,'r')
legend('Euler','R-K','Exact-soln')
end

function u = runge_kutta_44(f,u0,t,h,N)
u = zeros(N+1,1);
u(1) = u0;
for n = 1:N
    k1 = f(t(n),u(n));
    k2 = f(t(n)+0.5*h,u(n)+0.5*h*k1);
    k3 = f(t(n)+0.5*h,u(n)+0.5*h*k2);
    k4 = f(t(n)+h,u(n)+h*k3);
    u(n+1) = u(n) + h*( k1 + 2*k2 + 2*k3 + k4 )/6;
end
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 f = f1(t,u)
f = -5*u;
end

function f = exact1(t)
f = exp(-5*t);
end

function f = f2(t,u)
f = 4*t*sqrt(u);
end

function f = exact2(t)
f = ( 1 + t.^2 ).^2;
end