Code in Matlab: Complete lines s1,s2, s3, p(p=3) , yn1 and err so that the ode32
ID: 3834819 • Letter: C
Question
Code in Matlab: Complete lines s1,s2, s3, p(p=3) , yn1 and err so that the ode324 code implements the embedded explicit trapezium and simpson pair.
Code:
% INPUTS: A function handle f = @(t,y) ... such that dy/dt = f(t,y)
% The times span tspan(1) < t <tspane(end)
% The initial condition y(tspan(1)) - y0
% The required accuracy, tol.
% OUTPUTS: The selected time steps tt, and solution values yy.
% (Note, if no outputs requested, then a plot is shown instead).
% Initial step size
h = sqrt(tol);
% Initialise solver
t = tspan(1);
tend = tspan(2);
yn = y0;
flag = 0;
% Storage
tt = t;
yy = yn;
%% Loop through time
while ( t < tend )
% RK stages
s1 = ;
s2 = ;
s3 = ;
p = ;
yn1 = yn + ;
% Error estimate
err = ;
% Step size adjust:
hnew = 0.8*(tol*min(abs(yn1)./err))^(1/(p+1))*h;
hnew = min(hnew, 2*h); % Don't let h grow too fast.
if ( norm(err, inf) > tol ) % Error is too large:
if ( ~flag ) % Try the new time step,
flag = 1; % Set flag to remember failure.
else % Time step failed.
hnew = h/2; % Take half of old one,
end
else % Error is small enough!
% Store solution and move to next time:
yn = yn1;
t = t + h;
yy = [yy, yn1];
tt = [tt, t];
flag = 0; % Successful step. Set flag to happy.
end
% Choose next step (ensure we finish at the end of the interval!)
h = min(hnew, tend - t);
end
if ( nargout == 0 )
subplot(2,1,1)
plot(tt, yy(1,:), '.'), shg
title('First component of solution')
subplot(2,1,2)
semilogy(tt(2:end), abs(diff(tt)))
title(['No. of time steps = ' int2str(length(tt))])
tt = [];
yy = [];
end
end
Explanation / Answer
Code:
% INPUTS: A function handle f = @(t,y) ... such that dy/dt = f(t,y)
% The times span tspan(1) < t <tspane(end)
% The initial condition y(tspan(1)) - y0
% The required accuracy, tol.
% OUTPUTS: The selected time steps tt, and solution values yy.
% (Note, if no outputs requested, then a plot is shown instead).
% Initial step size
h = sqrt(tol);
% Initialise solver
t = tspan(1);
tend = tspan(2);
yn = y0;
flag = 0;
% Storage
tt = t;
yy = yn;
%% Loop through time
while ( t < tend )
% RK stages
s1 = ;
s2 = ;
s3 = ;
p = ;
yn1 = yn + ;
% Error estimate
err = ;
% Step size adjust:
hnew = 0.8*(tol*min(abs(yn1)./err))^(1/(p+1))*h;
hnew = min(hnew, 2*h); % Don't let h grow too fast.
if ( norm(err, inf) > tol ) % Error is too large:
if ( ~flag ) % Try the new time step,
flag = 1; % Set flag to remember failure.
else % Time step failed.
hnew = h/2; % Take half of old one,
end
else % Error is small enough!
% Store solution and move to next time:
yn = yn1;
t = t + h;
yy = [yy, yn1];
tt = [tt, t];
flag = 0; % Successful step. Set flag to happy.
end
% Choose next step (ensure we finish at the end of the interval!)
h = min(hnew, tend - t);
end
if ( nargout == 0 )
subplot(2,1,1)
plot(tt, yy(1,:), '.'), shg
title('First component of solution')
subplot(2,1,2)
semilogy(tt(2:end), abs(diff(tt)))
title(['No. of time steps = ' int2str(length(tt))])
tt = [];
yy = [];
end
end