Parallal computation: Splitting a huge ode45 function by building own runge-kutta-methods

Hi guys,

i hope you have a good idea. ;) I have to compute an imense ode45 function, with approx. 100 function parts and 50,000 time steps. Unfortenately, this takes on my computer several hours or days, so I tried some ideas to increase the speed. However, I do not want to use the dedicated matlab servers. To be presice, by '100 function parts', I mean this:

function dy = myODE(t,x)
dy = zeros(100,1)
...

Since each function step is dependent on all function values of the step before, the only way i see to use parallel programming is to split the function parts, and calculate only one time step before merging the results. Unfortunately, ode45 gives at least 42 function values and sub-steps, but I only need too. It is possible, to use only the very first and very last one, but this might be a waste of resources. Is there a way, to decrease the number of caluculated time steps?

However, I had the idea, to build my own range-kutta-method, based on the core skript parts of the original ode45 function. It looks like this:

function yOut = rungeKutta4(odeFcn,h,init)

B = [
    1/5         3/40    44/45   19372/6561      9017/3168       35/384
    0           9/40    -56/15  -25360/2187     -355/33         0
    0           0       32/9    64448/6561      46732/5247      500/1113
    0           0       0       -212/729        49/176          125/192
    0           0       0       0               -5103/18656     -2187/6784
    0           0       0       0               0               11/84
    0           0       0       0               0               0
    ];

E = [71/57600; 0; -71/16695; 71/1920; -17253/339200; 22/525; -1/40];

f = zeros(numel(init),7);
f(:,1) = init;

      hB = h*B;
      f(:,2) = feval(odeFcn,init+f*hB(:,1));
      f(:,3) = feval(odeFcn,init+f*hB(:,2));
      f(:,4) = feval(odeFcn,init+f*hB(:,3));
      f(:,5) = feval(odeFcn,init+f*hB(:,4));
      f(:,6) = feval(odeFcn,init+f*hB(:,5));
      ynew = init+f*hB(:,6);
      yOut = ynew';
      f(:,7) = feval(odeFcn,ynew);

%%Estimate the error.
err = h*(norm(f * E) / norm(ynew));
if err >= 1e-3
  disp('Error. Relative Tolerance greater than 1e-3')
end

end

[TheKawasakifan]

There are some changes on the code compared to ode45. First of all, my odeFunction is not time dependent, so I left that part completely out. Second, since only one time step is calculated, only the new values are returned as result. Third, the odeArguments part, the events part, the check error part and many more are left out, since i don't need them.

Of course, it did not work like it should. ;)

Here is the result of a simple swinging function:

blue is my range-kutta-method, red and correct is ode45.

The function which was evaluated was:

y' ' + y = 0

function testSwingParallel

init = [1; 0];
tstart = 0;
tend = 20;
h = 0.01;
nt = (tend - tstart)/h;
tspan = linspace(tstart,tend,nt);
solution = zeros(nt,numel(init));

initTemp = init;
tic
for k = 1:nt

      solution(k,:) = rungeKutta4(@myFcn,h,initTemp);
      initTemp = solution(k,:)';
  end
  timeRungeKutta = toc;

tic
[T,Y] = ode45(@myFcnOde,tspan,init);
timeOde45 = toc;

tspan = linspace(tstart,tend,nt);
plot(tspan,solution(:,1),'b')
hold on
plot(T,Y(:,1),'r')
hold off

disp(['timeRungeKutta = ', num2str(timeRungeKutta)]);
disp(['timeOde45 = ', num2str(timeOde45)]);

end

function dy = myFcn(x)

dy = zeros(2,1);
dy(2) = x(1);
dy(1) = -x(2);
end

function dy = myFcnOde(~,x)

dy = zeros(2,1);
dy(2) = x(1);
dy(1) = -x(2);
end

[TheKawasakifan]

I have also tried a simpler runge-kutta-code, it produced the same worng results.

function solution = rungeKutta4-2(fun,dt,x0)
ndeq = length (x0);
x = zeros (1,ndeq);
x(1,:) = x0';

      k1 = feval(fun,x(1,:));
      k2 = feval(fun,x(1,:)+k1'*dt/2);
      k3 = feval(fun,x(1,:)+k2'*dt/2);
      k4 = feval(fun,x(1,:)+k3'*dt);
      solution = x (1,:) + (k1' + 2*k2' + 2*k3' + k4') * dt/6 ;
  end

Any ideas have to solve this problem?