Fourth Order Runge-Kutta Method for the System of three Differential Equation

i want to solve that system using 4 order kutta method

i wrote that code but results are very weird, Q and R are two functions i wrote seperately , can someone see what's missing ?

clc
clear all
close all
%Constantes 
Mp=0.1;
Ap=(5.026548246)*10^(-3);
Patm=101325;
Kp=78750;
Xpo=877/78750;
rou=880;
Cd=0.7;
Aorifice=1.963495408*(10^(-5));
Vo=7*10^(-5);
Beta=1600*(10^5);
% Initialise step-size variables
h = 0.00001;
%t = (0:h:1)';
%N = length(t);
%x1=zeros(1,length(t));
%x2=zeros(1,length(t));
%x3=zeros(1,length(t));
% Conditions initiales
u=0;
x1(1)=0;
x2(1)=0;
x3(1)=0;
t(1)=0;
for i=1:100000
    K1P=x2(i);
    K1Q=Q(x1(i),x3(i));
    K1R=R(u(i),x1(i),x2(i),x3(i));
    
    K2P=x2(i)+(h/2)*K1Q;
    K2Q=Q(x1(i)+(h/2)*K1P,x3(i)+(h/2)*K1R);
    K2R=R(u(i),x1(i)+(h/2)*K1P,x2(i)+(h/2)*K1Q,x3(i)+(h/2)*K1R);
    
    
    K3P=x2(i)+(h/2)*K2Q;
    K3Q=Q(x1(i)+(h/2)*K2P,x3(i)+(h/2)*K2R);
    K3R=R(u(i),x1(i)+(h/2)*K2P,x2(i)+(h/2)*K2Q,x3(i)+(h/2)*K2R);
    
    
    K4P=x2(i)+h*K3Q;
    K4Q=Q(x1(i)+h*K3P,x3(i)+h*K3R);
    K4R=R(u(i),x1(i)+h*K3P,x2(i)+h*K3Q,x3(i)+h*K3R);
    
   
    x1(i+1)=x1(i)+(h/6)*(K1P+(2*K2P)+(2*K3P)+K4P);
    x2(i+1)=x2(i)+(h/6)*(K1Q+(2*K2Q)+(2*K3Q)+K4Q);
    x3(i+1)=x3(i)+(h/6)*(K1R+(2*K2R)+(2*K3R)+K4R);
    u(i+1)=2.995134969444502*(10^3)*(t(i)^10)+(2.454194242827056*(10^3))*(t(i)^9)-(3.074449680687996*(10^4))*(t(i)^8)-(1.831814183200164*(10^4))*(t(i)^7)+(1.225116448716314*(10^5))*(t(i)^6)+(1.225116448716314*(10^5))*(t(i)^5)-(2.265044070414277*(10^5))*(t(i)^4)+(6.259310423328208*(10^2))*(t(i)^3)+(1.720638633336634*(10^5))*(t(i)^2)+(8.587517861009773*(10^4))*t(i)+(2.102731998692876*(10^5));
    t(i+1)=t(i)+h;
    
end

subplot(2,2,1)
plot(t,x1)
xlabel('time')
ylabel('x1')


subplot(2,2,2)
plot(t,x2)
xlabel('time')
ylabel('x2')

subplot(2,2,3)
plot(t,x3)
xlabel('time')
ylabel('x3')

[Ameer Hamza]

Do you want to write your own code for RK4, or can you use MATLAB's function such as ode45.

[JEDIDI Ran]

i want to use my own code

Answer 1

[James Tursa]

I don't really want to sift through all of that code, mostly because you are using different variables for the various states. I would advise you vectorize everything into one state vector and write your derivative function based on that. This will cut down on your hand written code to 1/3 of what you have, make it simpler, and put it in a form that you can check results directly against ode45( ). E.g., define a 3-element state vector x and your individual states as:

x(1) = x1

x(2) = x2

x(3) = x3

Then write a derivative function for x, where x is a 3-element state vector:

function dxdt = myderiv(t,x,Mp,Ap,Kp,...etc) % you fill in all of the constants needed here
x1 = x(1);
x2 = x(2);
x3 = x(3);
u = _____;  % you fill this in
dx1dt = _____; % you fill this in
dx2dt = _____; % you fill this in
dx3dt = _____; % you fill this in
dxdt = [dx1dt;dx2dt;dx3dt];
end

Then in your RK4 code you would define the function handle:

f = @(t,x) myderiv(t,x,Mp,Ap,Kp,...etc) % again, you fill this in

And your RK4 code then becomes simplified since you only have one set of vectorized equations to write using that single f function handle. The time series of state vectors will be a matrix, with each column of the matrix being the state vector at a particular time. E.g., you would initialize your state as

x(:,1) = _____; % you fill this in, a 3x1 column vector for initial x1, x2, x3

and your looping RK4 code would have something like this:

    K1 = f( t(i)      , x(:,i)          );
    K2 = f( t(i) + h/2, x(:,i) + K1*h/2 );
    K3 = f( t(i) + h/2, x(:,i) + K2*h/2 );
    K4 = f( t(i) + h  , x(:,i) + K3*h   );
    x(:,i+1) = x(:,i) + (h/6)*( K1 + 2*K2 + 2*K3 + K4 );