lab12ex1.m

来自「matlab codes for solving ODEs, including」· M 代码 · 共 26 行

function lab12ex1
t0 = 0;
tf = 1;
y0 = 1;
yexact = exp(-1);
for h = [0.1,0.01]
    N = round((tf-t0)/h);
    [t,ye]=euler(@f,[t0,tf],y0,N);
    [t,yt]=taylor(@f,@ft,@fy,[t0,tf],y0,N);
    [t,yr]=rk4(@f,[t0,tf],y0,N);
    tt = t0:h:tf;
    plot(tt,ye,'r',tt,yt,'b',tt,yr,'k');
    disp(['with h=' num2str(h) ...
            ' Euler error=' num2str(yexact-ye(end)) ...
            'Taylor error=' num2str(yexact-yt(end)) ...
            '  RK-4 error=' num2str(yexact-yr(end))])
end

function dydt = f(t,y)
dydt = -10*(y-exp(-t)) - exp(-t);

function dfdt = ft(t,y)
dfdt = -9*exp(-t);

function dfdy = fy(t,y)
dfdy = -10;