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$
2
  • 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, 2016 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, 2016 at 0:36

1 Answer 1

Reset to default
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$
4
  • $\begingroup$ It does not work... I get another bad solution $\endgroup$
    – Matheus Danella
    Nov 2, 2016 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, 2016 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, 2016 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, 2016 at 0:55

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .