% stability lecture, June 10th, 2015 close, clear % Curtis-Hirschfelder example R = 50; R = 2000; f = @(t,y) -R*(y-cos(t)); N = 40; % number of timesteps %N = 80; %N = 120; h = 1.5/N; % step size t(1) = 0; % initial time % initial values w(1) = 0; y(1) = 0; z(1) = 0; for n=1:N-1 t(n+1) = t(n) + h; % explicit Euler %w(n+1) = w(n) + h*f(t(n),w(n)); % implicit Euler y(n+1) = fsolve(@(x) x-y(n)-h*f(t(n+1),x),y(n)); % trapeziodal rule z(n+1) = fsolve(@(x) x-z(n)-h/2*(f(t(n),z(n)) + f(t(n+1),x)),z(n)); end %hold on %plot(t,w,'r-d',t,y,'b-o',t,z,'g-*') plot(t,y,'b-o',t,z,'g-*') %legend('explicit Euler','implicit Euler', 'TR') legend('implicit Euler', 'TR')