Runge-Kutta function with a second order ODE

Hello, I have this section of code I am trying to work with. I need it to be able to accept a matrix form of a second order derrivative. I'll post an example below.

function [derriv_value] = FunctionC(x,y)
%Function that contains the derrivative value
derriv_value = [y(2); -9*y(1)+sin(x)]; % y(1) = y , y(2) = v
end

This is my function I am calling into my Runge-Kutta function. It is a second order ODE. I need my Runge-Kutta to be able to accept it, but I am not sure how. I tried altering how the inputs to the equation are formatted but nothing has worked. Here is the Runge-Kutta code.

function [x, yvecb] = MyVec_Function2(F,h,x0,x1,y0,y1)
% Note that F function expression is defined via Function Handle: F = @(x, y)(x+y)               
% change the function as you desire
                                             % step size (smaller step size gives more accurate solutions)
x = x0:h:x1;                                          % x space
yvecb(:,1)= y0;                         % Memory allocation
y(y0) = y1;                                       % initial condition
for i=1:(length(x)-1)                             
% i=1:(length(x)-1)                              % calculation loop
    k1 = F(x(i),y(:,i));
    k2 = F(x(i)+0.5*h,y(i)+0.5*h*k1);
    k3 = F((x(i)+0.5*h),(y(:,i)+0.5*h*k2));
    k4 = F((x(i)+h),(y(:,i)+k3*h));
    y(:,i+1) = y(:,i) + (1/6)*(k1+2*k2+2*k3+k4)*h;  % main equation
end
figure, plot(x, y)  % To see the solution results
end

[James Tursa]

First, the function handle would simply be defined as:

F = @FunctionC

But, more importantly, this looks like a boundary value problem. I.e, it looks like you have two values of y at two different x's as boundary conditions. Is this true? If so, RK4 is not the method for you since it is an initial value solver, not a boundary value solver.

The code you have written will not work even for an initial value problem because the initial conditions are incorrect, in particular this:

yvecb(:,1)= y0;
y(y0) = y1; 

[OvercookedRamen]

The inital values used to be correct. The un altered code looks like this. I am simply trying to use a second order ODE and enter it in the form of a matrix/vector rather than a 1 dimensional value. Does that help any?

function [x, y] = FunctionBeta_Executor(F,h,x0,x1,y0,y1)
% Note that F function expression is defined via Function Handle: F = @(x, y)(x+y)               
% change the function as you desire
                                             % step size (smaller step size gives more accurate solutions)
x = x0:h:x1;                                          % x space
y = zeros(1,length(x));                             % Memory allocation
y(y0) = y1;                                       % initial condition
for i=1:(length(x)-1)                             
% i=1:(length(x)-1)                              % calculation loop
    k1 = F(x(i),y(i));
    k2 = F(x(i)+0.5*h,y(i)+0.5*h*k1);
    k3 = F((x(i)+0.5*h),(y(i)+0.5*h*k2));
    k4 = F((x(i)+h),(y(i)+k3*h));
    y(i+1) = y(i) + (1/6)*(k1+2*k2+2*k3+k4)*h;  % main equation
end
figure, plot(x, y)  % To see the solution results
end
% To run this function, supply the needed values. 
% F, can be any function file, example @FunctionA.
% Example, FunctionBeta_Executor(@FunctionA.....)
% h is the step size, x0 is the first x value, x1 is the second
% y0 and y1 supply the inital conditions, example, y(0)=1, y0 = 0, y1 = 1.

Answer 1

[James Tursa]

I'm still confused about your initial conditions and what y0 and y1 are. But let's back up a bit.

For a 2nd order ODE, your state vector will have 2 elements. E.g., if y is your state then y(1) could be position and y(2) could be velocity. This looks like the setup you have in your FunctionC( ) above. So for that case, the states in your RK4 solver have to reflect this and be 2-element vectors, not scalars. So that y pre-allocation needs to have 2 rows, not 1. And all the downstream indexing needs to account for this as well. E.g., something like this if we assume that those y0 and y1 values are the initial values at the x0 time (CAUTION: Not checked):

function [x, y] = FunctionBeta_Executor(F,h,x0,x1,y0,y1)
% Note that F function expression is defined via Function Handle: F = @(x, y)(x+y)               
% change the function as you desire
                                             % step size (smaller step size gives more accurate solutions)
x = x0:h:x1;                                          % x space
y = zeros(2,length(x)); % CHANGED                            % Memory allocation
y(1,1) = y0;  % CHANGED                                       % initial condition
y(2,1) = y1;  % NEW
for i=1:(length(x)-1)                             
% i=1:(length(x)-1)                              % calculation loop
    k1 = F( x(i)         , y(:,i)            ); % CHANGED
    k2 = F( x(i) + 0.5*h , y(:,i) + 0.5*h*k1 ); % CHANGED
    k3 = F( x(i) + 0.5*h , y(:,i) + 0.5*h*k2 ); % CHANGED
    k4 = F( x(i) +     h , y(:,i) +     h*k3 ); % CHANGED
    y(:,i+1) = y(:,i) + (1/6)*( k1 + 2*k2 + 2*k3 + k4 )*h;  % main equation % CHANGED
end
figure, plot(x, y(1,:))  % To see the solution results % CHANGED
end

Note also how much easier it is to read (and debug) when you use spacing!

[OvercookedRamen]

Yeah that seemed to do it. Thanks for explaing y(1) and y(2). I get this stuff usually pretty easy when I do by hand, but when it comes to putting it into MatLab it just confuses me so much. I'll do my best to keep my code clean for now on. Thanks a lot man. I ran it, it works like a charm.

Answer 2

[khalida]

clear,clc

% % % % % RK4 system of odes higher order

h = 0.1;

x_beg = 0;

x_end = 0.7;

y_initial= 0;

z_initial= 0;

F_xy=@(x,z)(z);

F_xz=@(x,y,z)(-((1+5*x)/2*x*(x+1))*z-y^3+5*x+11*x^2+(8/27)*x^9);

x=x_beg:h:x_end;

y=zeros(1,length(x));

y(1)=y_initial;

z=zeros(1,length(x));

z(1)=z_initial;

for i=1:(length(x)-1);

i=1:(length(x)-1);

k1 = F_xy(x(i),y(i),z(i));

L1 = F_xz(x(i),y(i),z(i));

k2 = F_xy((x(i)+h/2),(y(i)+0.5*h*k1),z(i)+0.5*h*L1);

L2 = F_xz((x(i)+h/2),(y(i)+0.5*h*k1),z(i)+0.5*h*L1);

k3 = F_xy((x(i)+h/2),(y(i)+0.5*h*k2),z(i)+0.5*h*k2);

L3 = F_xz(x(i)+h/2,(y(i)+0.5*h*k2),(z(i)+0.5*h*L2));

k4 = F_xy((x(i)+h),(y(i)+k3*h),z(i)+L3*h);

L4 = F_xz((x(i)+h),(y(i)+k3*h),(z(i)+L3*h));

y(i+1) = y(i) + (1/6)*(k1+2*k2+2*k3+k4)*h;

z(i+1) = z(i) + (1/6)*(L1+2*L2+2*L3+L4)*h;

end

% plot(y,z)

figure(1)

plot(x,z)

hold on

plot(x,y)

Answer 3

[khalida]

clear,clc

% % % % % RK4 system of odes higher order

h = 0.1;

x_beg = 0;

x_end = 0.7;

y_initial= 0;

z_initial= 0;

F_xy=@(x,z)(z);

F_xz=@(x,y,z)(-((1+5*x)/2*x*(x+1))*z-y^3+5*x+11*x^2+(8/27)*x^9);

x=x_beg:h:x_end;

y=zeros(1,length(x));

y(1)=y_initial;

z=zeros(1,length(x));

z(1)=z_initial;

for i=1:(length(x)-1);

i=1:(length(x)-1);

k1 = F_xy(x(i),y(i),z(i));

L1 = F_xz(x(i),y(i),z(i));

k2 = F_xy((x(i)+h/2),(y(i)+0.5*h*k1),z(i)+0.5*h*L1);

L2 = F_xz((x(i)+h/2),(y(i)+0.5*h*k1),z(i)+0.5*h*L1);

k3 = F_xy((x(i)+h/2),(y(i)+0.5*h*k2),z(i)+0.5*h*k2);

L3 = F_xz(x(i)+h/2,(y(i)+0.5*h*k2),(z(i)+0.5*h*L2));

k4 = F_xy((x(i)+h),(y(i)+k3*h),z(i)+L3*h);

L4 = F_xz((x(i)+h),(y(i)+k3*h),(z(i)+L3*h));

y(i+1) = y(i) + (1/6)*(k1+2*k2+2*k3+k4)*h;

z(i+1) = z(i) + (1/6)*(L1+2*L2+2*L3+L4)*h;

end

% plot(y,z)

figure(1)

plot(x,z)

hold on

plot(x,y)