Borrar filtros
Borrar filtros

Runge-Kutta function with a second order ODE

104 visualizaciones (últimos 30 días)
OvercookedRamen
OvercookedRamen el 29 de Oct. de 2019
Respondida: khalida el 17 de Jul. de 2023
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 comentarios
James Tursa
James Tursa el 29 de Oct. de 2019
Editada: James Tursa el 29 de Oct. de 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 el 29 de Oct. de 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.

Iniciar sesión para comentar.

Respuesta aceptada

James Tursa
James Tursa el 29 de Oct. de 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 comentario
OvercookedRamen
OvercookedRamen el 29 de Oct. de 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.

Iniciar sesión para comentar.

Más respuestas (2)

khalida
khalida el 14 de Jul. de 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 el 17 de Jul. de 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)

Categorías

Más información sobre Programming en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by