the stiff IVP y'(t) = -30*y+30*t^2+2*t; y(0)=1, & exact solution is y(t)=e^-30*x +x^2. using fourth order Runge-kutta method. solve by hand and by matlab

function rungekutta h = 0.5; t = 0; w = 0.5; fprintf('Step 0: t = %12.8f, w = %12.8f\n', t, w); for i=1:4 k1 = h*f(t,w); k2 = h*f(t+h/2, w+k1/2); k3 = h*f(t+h/2, w+k2/2); k4 = h*f(t+h, w+k3); w = w + (k1+2*k2+2*k3+k4)/6; t = t + h; fprintf('Step %d: t = %6.4f, w = %18.15f\n', i, t, w); end y_exact=exp(-30.*t)+t.^2; subplot(2,1,1) plot(t,w-y_exact) ylabel('exact vs. approximated function') subplot(2,1,2) plot(t,w-y_exact) ylabel('error') %%%%%%%%%%%%%%%%%% function v = f(t,y) v = -30*y+30*t+2*t;

Answer 1

[Mischa Kim]

Jamila, please format your question(s) for better readability.

There are a couple of issues with your code. Check out the following:

function rungekutta 

h     = 0.01; 
tend  = 5;
tspan = 0:h:tend;
t(1)  = tspan(1); 
w(1)  = 0.5; 

fprintf('Step 0: t = %12.8f, w = %12.8f\n', t, w); 
for ii=1:length(tspan) 
    k1 = h*f(t(ii),w(ii)); 
    k2 = h*f(t(ii)+h/2, w(ii)+k1/2);
    k3 = h*f(t(ii)+h/2, w(ii)+k2/2); 
    k4 = h*f(t(ii)+h, w(ii)+k3); 
    w(ii+1)  = w(ii) + (k1+2*k2+2*k3+k4)/6; 
    t(ii+1)  = t(ii) + h; 
    %fprintf('Step %d: t = %6.4f, w = %18.15f\n', ii, t(ii), w(ii)); 
end
y_exact=exp(-30.*tspan)+tspan.^2; 
subplot(2,1,1) 
plot(tspan,w(1:length(w)-1),tspan,y_exact) 
ylabel('exact vs. approximated function') 
subplot(2,1,2)
plot(tspan,w(1:length(w)-1)-y_exact) 
ylabel('error') %%%%%%%%%%%%%%%%%%function v = f(t,y) v = -30*y+30*t+2*t;
end

function f_eval = f(t,y)
  f_eval = -30*y + 30*t^2 + 2*t;
end

• In the loop you need to save each time step data (t,w) as matrices so you can access (plot) the graph later on.
• A time step of h = 0.5 is huge, especially for an unbounded function.
• For f(t(ii),w(ii)) you need to define the function accordingly.