0
$\begingroup$

I have been trying to solve a system of 8 coupled ODEs using the 4th order Runge-Kutta method. Also, I am using Matlab to write the algorithm. The ODE system can be seen below:

System of 8 coupled ODEs

The initial conditions are:

Initial Conditions

Also, V1 dot and V1 double dot are zero since. Furthermore, beta is just a constant parameter, since I want to make the code work it can be choosen as 1.

I have been checking other topics here about the Runge-Kutta method, and I have got some nice ideas about how to apply it on Matlab, and I wrote some chunk of code, it can also be seen below:

% Parâmeters
Beta = 2;
q = 0.1         % Iteration intervals
qfinal = 55;  

% Initial conditions
x1(1)=0;
x2(1)=0;
x3(1)=0;
x4(1)=0;
x5(1)=0;
x6(1)=0;
x7(1)=0;
x8(1)=0;


% Definir a EDO

fx1=@(x2) x2;
fx2=@(x1,x3,x5) (-7*Beta)*x2 + (4*Beta)*x3 + (-1*Beta)*x5;
fx3=@(x4) x4;
fx4=@(x1,x3,x5,x7) (4*Beta)*x1 + (-6*Beta)*x3 + (4*Beta)*x5 + (-1*Beta)*x7;;
fx5=@(x6) x6;
fx6=@(x1,x3,x5,x7) (-1*Beta)*x1 + (4*Beta)*x3 + (-5*Beta)*x5 + (2*Beta)*x7;;
fx7=@(x8) x8;
fx8=@(x3,x5,x7) (-2*Beta)*x3 + (4*Beta)*x5 + (-2*Beta)*x7;

t =1:ceil(qfinal/q)+1;
% RK4 Loop
for i=1:ceil(qfinal/q)

    x1_1=fx1(x2(i));
    x2_1=fx2(x1(i),x3(i),x5(i));
    x3_1=fx3(x4(i));
    x4_1=fx4(x1(i),x3(i),x5(i),x7(i));
    x5_1=fx5(x6(i));
    x6_1=fx6(x1(i),x3(i),x5(i),x7(i));
    x7_1=fx7(x8(i));
    x8_1=fx8(x3(i),x5(i),x7(i));

    x1_2=fx1(x2(i)+0.5*x2_1);
    x2_2=fx2(x1(i)+0.5*q,x3(i)+0.5*q*x3_1,x5(i)+0.5*q*x5_1);
    x3_2=fx3(x4(i)+0.5*q);
    x4_2=fx4(x1(i)+0.5*q,x3(i)+0.5*q*x3_1,x5(i)+0.5*q*x5_1,x7(i)+0.5*q*x7_1);
    x5_2=fx5(x6(i)+0.5*q);
    x6_2=fx6(x1(i)+0.5*q,x3(i)+0.5*q*x3_1,x5(i)+0.5*q*x5_1,x7(i)+0.5*q*x7_1);
    x7_2=fx7(x8(i)+0.5*q);
    x8_2=fx8(x3(i)+0.5*q,x5(i)+0.5*q*x5_1,x7(i)+0.5*q*x7_1);

    x1_3=fx1(x2(i)+0.5*q);
    x2_3=fx2((x1(i)+0.5*q),(x3(i)+0.5*q*x3_2),(x5(i)+0.5*q*x5_2));
    x3_3=fx3(x4(i)+0.5*q);
    x4_3=fx4((x1(i)+0.5*q),(x3(i)+0.5*q*x3_2),(x5(i)+0.5*q*x5_2),(x7(i)+0.5*q*x7_2));
    x5_3=fx5(x6(i)+0.5*q);
    x6_3=fx6((x1(i)+0.5*q),(x3(i)+0.5*q*x3_2),(x5(i)+0.5*q*x5_2),(x7(i)+0.5*q*x7_2));
    x7_3=fx7(x8(i)+0.5*q);
    x8_3=fx8((x3(i)+0.5*q),(x5(i)+0.5*q*x5_2),(x7(i)+0.5*q*x7_2));

    x1_4=fx1(x2(i)+q);
    x2_4=fx2((x1(i)+q),(x3(i)+q*x3_3),(x5(i)+q*x5_3));
    x3_4=fx3(x4(i)+q);
    x4_4=fx4((x1(i)+q),(x3(i)+q*x3_3),(x5(i)+q*x5_3),(x7(i)+q*x7_3));
    x5_4=fx5(x6(i)+q);
    x6_4=fx6((x1(i)+q),(x3(i)+q*x3_3),(x5(i)+q*x5_3),(x7(i)+q*x7_3));
    x7_4=fx7(x8(i)+q);
    x8_4=fx8((x3(i)+q),(x5(i)+q*x5_3),(x7(i)+q*x7_3));

    x1(i+1) = x1(i) + (1/6)*(x1_1+2*x1_2+2*x1_3+x1_4)*q;
    x2(i+1) = x2(i) + (1/6)*(x2_1+2*x2_2+2*x2_3+x2_4)*q;
    x3(i+1) = x3(i) + (1/6)*(x3_1+2*x3_2+2*x3_3+x3_4)*q;
    x4(i+1) = x4(i) + (1/6)*(x4_1+2*x4_2+2*x4_3+x4_4)*q;
    x5(i+1) = x5(i) + (1/6)*(x5_1+2*x5_2+2*x5_3+x5_4)*q;
    x6(i+1) = x6(i) + (1/6)*(x6_1+2*x6_2+2*x6_3+x6_4)*q;
    x7(i+1) = x7(i) + (1/6)*(x7_1+2*x7_2+2*x7_3+x7_4)*q;
    x8(i+1) = x8(i) + (1/6)*(x8_1+2*x8_2+2*x8_3+x8_4)*q;
end

Then, the code plots the solutions

figure
plot(t,x1,t,x3,t,x5,t,x7)

figure
plot(t,x2,t,x4,t,x6,t,x8)

The problem is that I should get oscillating functions as solutions, but I am getting exponential ones. I am pretty sure something is wrong in my code, but I am not sure of what it is.

A mistaken velocity definition on Matlab might the be cause, but I am not sure about anything, and that is why I am asking for help here.

$\endgroup$
  • 1
    $\begingroup$ You should investigate the use of vectors, vector-valued functions and the utility of matrix-vector products when using the MATrix LABoratory. $\endgroup$ – Lutz Lehmann Nov 2 '16 at 16:59
  • $\begingroup$ As you have a linear system $X'(t)=AX(t)$, with $A$ a $8 \times 8$ matrix, why don't you directly use the explicit solution $X(t)=expm(tA)*X(0)$? (expm = matrix exponential). It works very well with Matlab. No need in this case to use Runge Kutta... $\endgroup$ – Jean Marie Nov 3 '16 at 0:36
1
$\begingroup$

Should

x8_4=fx8((x3(i)+q),(x5(i)+q*x5_3),(x7(i)+q*x7_3));

have an additional multiplicand like

x8_4=fx8((x3(i)+q*x3_3),(x5(i)+q*x5_3),(x7(i)+q*x7_3));

and similarly for other first-arguments throughout the code?

$\endgroup$
  • $\begingroup$ It does not work... I get another bad solution $\endgroup$ – Matheus Danella Nov 2 '16 at 18:14
  • $\begingroup$ You should get no non-trivial solution at all, since you start in the fixed point, the zero vector. The whole solution should be zero. -- Since you do not get that, you have not corrected the above error in all places from x2_2=fx2(x1(i)+0.5*q*x1_1,... on. $\endgroup$ – Lutz Lehmann Nov 2 '16 at 18:45
  • $\begingroup$ Yes I did!! May I send you the corrected code somehow? I mean, I did and I get a solution of zeroes only. $\endgroup$ – Matheus Danella Nov 7 '16 at 0:43
  • $\begingroup$ I think I get these because of the initial definitions, since they all say that x1_1...x8_1 are zero in the beginning. $\endgroup$ – Matheus Danella Nov 7 '16 at 0:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.