My code uses the Euler scheme to solve an equation.How would I convert it to Runge Kutta scheme to make it faster for convergence.

My code tries to solve the Cummins Equation using an Impluse Response Function in time domain. For a given external forcing, the resultant displacement against time is plotted. The equation seems to be stiff and is too slow and doesn't seem to converge.How do I convert it to Runge Kutta?

[John D'Errico]

Runge-Kutta may not be the very best choice for stiff systems either. Why not use the ODE solvers already in MATLAB? There are stiff solvers designed to solve that class of problem.

Recreating the wheel is just a bad idea in software. If high quality code already exists, then use it.

[Vignesh Kumar]

Thanks John. Thanks for your reply, I do not know how to use ode45 directly to solve this as my equation consists of an unknown velocity term.

Answer 1

[Ced]

Hi

Not sure if I understand the question correctly, but you basically already have the answer, and most of the code you need.

Assuming you want to implement your own version of Runge Kutta and not use the matlab ode, you just need to slightly change 'myfun':

Yh=[z0;zd0];
for i=2:length(time)
   dY = myode(time(i),dt,Yh);
   temp=Yh(:,i-1)+dY*dt; %Euler integration
   Yh=cat(2,Yh,temp);
end

Just replace the line with % Euler integration by the Runge Kutta scheme of your choice. The most popular one is RK4, you'll find everything you need on wikipedia: Runge Kutta Methods

One suggestion though: Instead of concatenating your result as Yh=cat(2,Yh,temp), you should preallocate the whole vector, e.g.:

Y = zeros(2,length(time));
Y(:,1) = [z0;zd0];
for i = 2:length(time)
   ...
   Yh(:,i) = Yh(:,i-1) + ...; <-- plug in update here, i.e. Yh(:,i-1) + dt*(k1+2*k2+2*k3+k4)/6 for RK4
end

[Vignesh Kumar]

Thanks Ced.Thanks for your reply. My code works with Euler, but I do not see any convergence in the solution. Irrespective of the forcing I use, I get similar results. And the result is inconsistent as I increase time span. SO I wonder if my code is right? or my code is fine and the numerical scheme is not robust.? I would be extremely grateful if you could help me with modifying the code to achieve a solution. My doubt is that after i remove the impulse momentary force, the result should gradually decay to zero, which doesn't happen in my code.

[Ced]

You can just select a very small timestep for your Euler scheme. If things change, then the integration could be a problem. It looks like you have some discontinuities in your equations though, so the choice of discretization can have quite an impact, even if things are working correctly.

Your can also try implementing the RK4 scheme as you suggested. (Just replace the Euler integration line by the RK4 scheme you will find on the wiki page). If the result is similar as with Euler, then you know it's your code and not the integration.

Answer 2

[Vignesh Kumar]

I have replaced the code solver as RK4 and find some better results, showing some hope of convergence.Very small time step or very large time takes a lot of run time still, but the solution looks much better by RK than Euler, So is my code right? Ive attached both results by Euler and RK4 as screenshots.