1

I have a Matlab function for doing Runge-Kutta4k approximation for first-order ODE's, and I want to adapt it to work for second-order ODE's. Would anyone be able to help me get started? This is what I have for first order RK4K:

function [y,t,h] = rungekutta4kh(F,y0,a,b,n)

% Euler ODE solver
t = linspace(a,b,n);
h = t(2)-t(1);

y = zeros(n,1);
y(1) = y0;
for i=2:n
    s1 = feval(F,t(i-1),y(i-1));
    s2 = feval(F,t(i-1)+h/2,y(i-1)+ (h/2)*s1);
    s3 = feval(F,t(i-1)+h/2,y(i-1)+ (h/2)*s2);
    s4 = feval(F,t(i-1)+h,y(i-1)+ h*s3);

    y(i) = y(i-1) + ...
        (h/6)*( s1 + 2*s2 + 2*s3 + s4 );
end
kathystehl
  • 831
  • 1
  • 9
  • 26

1 Answers1

0

You transform it into a first order system, which means that the function F becomes vector valued and you should include the shape of y0 in the construction of the y list.

y''=f(x,y,y')

gets transformed to

y0'=y1
y1'=f(x,y0,y1)

so that (with matlab indexing)

F(x,y) = [ y(2), f(x,y(1),y(2)) ]
Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • So, you're saying that I should transform the second-order ODE that I'm working with into a system of first-order ODE's and then apply the algorithm that I already have but set it up inputs for the 2-equation system of ODE's? – kathystehl Dec 14 '15 at 20:02
  • Yes, exactly that. There are methods for second order systems, Beeman made a comparative study of them for molecular dynamics, but they mostly apply to (energy) conservative systems. – Lutz Lehmann Dec 16 '15 at 13:59