Solving a system of ODE with Crank Nicholson

I have a set of ODEs as follows :-

d(A)/dt = f(A,B,V)

d(B)/dt g(A,B,V)

d(C)/dt = h(A,B,C,V)

d(D)/dt = i(A,B,C,V)

V = j(t)

Initial values are known and I need to solve the system using Crank-Nicholson Method. Pseudocode will be welcome.

[Alan Stevens]

The Crank-Nicolson method is a PDE solution technique. You seem to have a system of ODEs. Have a look at

doc ODE45

[Arkajyoti Chaterjee]

I know that but remain interested to solve this using CN (all ODEs can be treated as PDEs).

Also, if I do use ode45, how to pass 3/4 variable params? [Only 2 input arguments are accepted and this example passes constants.]

[Alan Stevens]

Read the documentation on ode45 to see how to pass more than 2 input arguments.

[Arkajyoti Chaterjee]

That's resolved, thanks! ODE45 is giving the correct solutions.

Answer 1

[Bjorn Gustavsson]

Note that you have a set of coupled ODEs not PDEs. The Crank-Nicholson method is developed for solving PDEs where we have both "temporal" and "spatial" derivatives, whereas you only seem to have tetmporal derivatives. This makes C-N a rather unusual choise, but is you see your ODEs as "single-position PDEs" you might simply use the equations in that wikipedia-article and scrap all spatial derivatives and boundary-condition-concerns. It will take a bit of time to write out the derivatives and the system of equations you'll need to solve. But it might be well worth for illustrating the method.

Otherwise just use the standard methods for solving coupled ODEs available in matlab (ode45, etc).

HTH

[Arkajyoti Chaterjee]

In case of ODE, the CN reduces to this (warning:pdf file) but am having trouble setting up the code for the implicit solving.

[Bjorn Gustavsson]

It is easier to use the odeNN functions of matlab instead of creating a specialized version of C-N to handle a simple set of ODEs. (If you didn't have right-hand sides that were such general functions of the solution it would be another matter). I suggest this:

function dABCDdt = ode4ABCD(t,ABCD,other,parameters)
V = j(t); % as per your case
A = ABCD(1);
B = ABCD(2);
C = ABCD(3);
D = ABCD(4);
% Change as necessary
dABCDdt = [f_fcn(A,B,C,D,V);g_fcn(A,B,C,D,V);h_fcn(A,B,C,D,V);i_fcn(A,B,C,D,V)];
end

This you then integrate with (for example) ode45:

ABCD0 = [1 2 3 4];  % Adjust
t_span = [0 10*pi]; % Adjust...
[t_out,ABCD] = ode45(@(t,abcd) ode4ABCD(t,abcd,1i,[atan(pi/4),exp(-1)]),t_span,ABCD0);

HTH

[Arkajyoti Chaterjee]

ODE45 is solved; thanks. The issue with CN is that I have been asked to use the particular algorithm.

[Bjorn Gustavsson]

Well that will give you a set of four non-linear equations - since you'll have to replace A, B, C and D on the right-hand sides with the average values of the values at the current time-step and the next time-step. That is:

should become something like:

From there you might get away with linearising the RHS or you might have to solve that nonlinear system of equations in A, B, C and D at i+1. Perhaps you might get a good enough understanding if you implement a 1-D problem first - that would give you some hands-on experience of where this goes programming-wise.

HTH