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.
5 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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?
2 comentarios
John D'Errico
el 25 de Mzo. de 2016
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.
Respuestas (2)
Ced
el 25 de Mzo. de 2016
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
2 comentarios
Ced
el 25 de Mzo. de 2016
Editada: Ced
el 25 de Mzo. de 2016
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.
Ver también
Categorías
Más información sobre Numerical Integration and Differential Equations en Help Center y File Exchange.
Productos
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!