function [] = prob_2_10_pt3 % EX_2_10_PT3 - % % Solve y' = f(t,y) with the trapezoidal rule: % % ynp1 = y_n + 0.5 * h * ( f(tn,yn) + f(tnp1,ynp1) ); % % using quasi-newtons method. % % Written by: % -- % John L. Weatherwax 2005-03-10 % % email: wax@alum.mit.edu % % Please send comments and especially bug reports to the % above email address. % %----- a = 0; b = 1; h = 0.1; nSteps = round((b-a)/h); ys = [ 0; 1]; for ii=1:nSteps, tn = a + (ii-1)*h; yn = ys(:,end); [tnp1,ynp1] = AM2ni(tn,yn,h,@ode,@dode); ys = [ ys, ynp1 ]; end ts = a + h*(0:nSteps); yTruth = [ sin(ts); cos(ts) ]; figure; subplot(1,2,1); plot( ts, ys(1,:), '-r', ts, yTruth(1,:), '-g' ); title( 'y1 approx & truth' ); subplot(1,2,2); plot( ts, ys(2,:), '-r', ts, yTruth(2,:), '-g' ); title( 'y2 approx & truth' ); return; %------------------------------------------------------------------------ %------------------------------------------------------------------------ function dydt = ode(t,y) dydt = [ 0 1; -1 0 ] * y; %------------------------------------------------------------------------ %------------------------------------------------------------------------ function J = dode(t,y) J = [ 0 1; -1 0 ]; %------------------------------------------------------------------------ %------------------------------------------------------------------------ function [tnp1,ynp1] = AM2ni(tn,yn,h,f,dfdy) % AM2NI - % % Written by: % -- % John L. Weatherwax 2005-03-10 % % email: wax@alum.mit.edu % % Please send comments and especially bug reports to the % above email address. % %----- % Since this is the quasi-newton method we evaluate the % derivative at yn only once and use it for all iterations: J = feval(dfdy,tn,yn); M = eye(size(J)) - 0.5*h*J; [L,U]=lu(M); % Initialize: tnp1 = tn+h; nIters = 0; ynp1 = yn; % Do one iteration outside the while loop: b = yn + 0.5 * h * feval(f,tn,yn) + 0.5 * h * feval(f,tnp1,ynp1) - ynp1; deltam = U \ ( L \ b ); ynp1 = ynp1 + deltam; % Do remaining iterations (if needed): nIters = 1; while( (norm(deltam) >= 1e-3 * norm(ynp1)) && nIters < 101 ) b = yn + 0.5 * h * feval(f,tn,yn) + 0.5 * h * feval(f,tnp1,ynp1) - ynp1; deltam = U \ ( L \ b ); ynp1 = ynp1 + deltam; nIters = nIters+1; end nIters return;