Solving nonlinear scalar ode's

I have to solve an ivp of the form: du/dt = f(t,u), u(t0) = u0, for a function u = u(t)?R for t=>t0, using Adams Bashforth 2nd order linear multistep method.

I have to write a function file with input t0 = initial time, tf = final time, u01=[u0; u1], n = number of time steps.

And output: t - is vector of length n+1 containing the times t0,t1,...,tn , and u - is a vector of length n+1 with the first element of u being u0 and with the (i+1) the element of u being the approximation ui for i=1,...,n.

So far i have: funtion [t,u]=ivpab2(t0,tf,u01,n) and i know that h=(tf-t0)/n and t(i)=t0+i*h The testing is to be done for the function u'=f(t,u)=-2*u+3*exp(-2*t)cos(3*t) which also needs to be implemented into this function file.

I would be extremely thankful for any help to how to even start with this.

Thank you.

Answer 1

[Matt Tearle]

Given that you have to use AB2 (with, apparently, a fixed stepsize) I going to take a guess that this is a numerical analysis homework problem? (If not, use one of the built-in MATLAB linear multistep methods, such as ode113.)

For fixed stepsize, you can calculate all the t values immediately. Look at either the linspace function or the range operator (:).

You can and should make room for the u values before you start to calculate them. Look at zeros.

Set the first value of u from the initial condition.

Loop to fill in the rest of the values, based on the previous ones. Given that you know how many steps you're going to take, use a for-loop.

You can hard-code the rate equation inside your loop if you like. But if you want to be a bit fancier, use an anonymous function handle to define it. Something like f = @(x,y) -2*y + ... Then you can evaluate it anywhere later by simply calling f(t(k),y(k)) (or whatever).

[Rebecca]

Thanks for your response Matt. I will give this problem a go (with your help) and let you know if I come across any more difficulties.

[Rebecca]

When implementing the Runge-Kutta order 4 scheme(rk4) u1,...,un are

computed using the algorithm:

for i = 0,1,...,n-1

k1 = f (ti, ui),

k2 = f (ti + h/2, ui + h/2*k1),

k3 = f (ti + h/2, ui + h/2*k2),

k4 = f (ti + h, ui + h*k3),

ui+1 = ui + h/6 (k1 + 2k2 + 2k3 + k4).

Would this be written in matlab as the following:

for i = 1:n

k1 = h * feval ( f, u(:,i) );

k2 = h * feval ( f, u(:,i) + k1/2 );

k3 = h * feval ( f, u(:,i) + k2/2 );

k4 = h * feval ( f, u(:,i) + k3 );

u(:,i+1) = u(:,i) + h*( k1 + 2*k2 + 2*k3 + k4 ) / 6;

end;

I'm not too sure if the "h" in the last line (u(:,i+1)) should be there are not.

Any help would be very much appreciated.