Runge-Kutta function with a second order ODE

241 views (last 30 days)
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
  2 Comments
James Tursa
James Tursa on 29 Oct 2019
Edited: James Tursa on 29 Oct 2019
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
OvercookedRamen on 29 Oct 2019
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.

Sign in to comment.

Accepted Answer

James Tursa
James Tursa on 29 Oct 2019
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!
  1 Comment
OvercookedRamen
OvercookedRamen on 29 Oct 2019
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.

Sign in to comment.

More Answers (2)

khalida
khalida on 14 Jul 2023
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)

khalida
khalida on 17 Jul 2023
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)

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!