Runge-Kutta function with a second order ODE
86 views (last 30 days)
Show older comments
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
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;
Accepted Answer
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!
More Answers (2)
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)
0 Comments
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)
0 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!