Matlab error in ode45 or fourth-order Runge-Kutta method to solve a system of coupled ODEs

Question

I am a beginner at Matlab programming and with the Runge-Kutta method as well.

I'm trying to solve a system of coupled ODEs using a 4th-order Runge-Kutta method for my project work.

here is my problem...

G = 1.4;
g = 1.4;
k = 0;
z = 0;
b = 0.166667;

syms n; 
x2 = symfun(sym('x2(n)'),[n]);
x1 = symfun(sym('x1(n)'),[n]);
x3 = symfun(sym('x3(n)'),[n]);
x4 = symfun(sym('x4(n)'),[n]);
x5 = symfun(sym('x5(n)'),[n]);

k1 = [x2 * x1 *n *(1 - z * x2)*(x1 - n) - 2 * x3 * n *(1 - z * x2) - x4^2 * x2 *(1 - z * x2)- G *x3 *x2 ]./ [( G * x3 - (x1 - n)^2 * x2 *(1 - z * x2)) * n];   
k2 = [x2 * (1 - z * x2)*(x1 * x2 * ( x1 - 2 *n)*( x1 - n) + 2* x3 * n + x4^2 * x2 ) ]./ [( G * x3 - (x1 - n)^2 * x2 *(1 - z * x2)) * n * (x1 - n)];
k3 = [x3 * x2 * (2 * n * x1 - n)^2 * ( 1 - z * x2) + G * x1 * (x1 - 2 *n)* (x1 - n) + x4^2 * G]./ [( G * x3 - (x1 - n)^2 * x2 *(1 - z * x2)) * n * (x1 - n)];
k4 = [x4 * ( x1 + n)] ./ [n * (x1- n)];
k5 = - [x5] ./ [n * (x1- n)];

f = @(n,x) [k1;  k2;  k3;  k4; k5];
[n,xa] = ode45(f,[0 1],[1-b 1/b 1-b 0.01 0.02]);

errors are

Error using odearguments (line 93)
@(N,X)[K1;K2;K3;K4;K5] returns a vector of length 1, but the length
of initial conditions vector is 5. The vector returned by @(N,X)[K1;K2;K3;K4;K5] and the initial conditions vector must have
the same number of elements.

Error in ode45 (line 114)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs,
odeFcn, ...

Please guide me how can I solve the above problem with fourth-order Runge-Kutta method...

Answer 1

The error stems from using symbolic functions (mixed with a function handle) with a numeric solver. You need to create numeric functions for ode45 to function properly (I also replaced all of the [ and ] with ( and ) for grouping):

G = 1.4;
g = 1.4;
k = 0;
z = 0;
b = 0.166667;

k1 = @(n,x) (x(2) * x(1)*n *(1 - z * x(2))*(x(1) - n) - 2 * x(3) * n *(1 - z * x(2)) - x(4)^2 * x(2) *(1 - z * x(2))- G *x(3) *x(2)) ./ (( G * x(3) - (x(1) - n)^2 * x(2) *(1 - z * x(2))) * n);   
k2 = @(n,x) (x(2) * (1 - z * x(2))*(x(1) * x(2) * ( x(1) - 2 *n)*( x(1) - n) + 2* x(3) * n + x(4)^2 * x(2) )) ./ (( G * x(3) - (x(1) - n)^2 * x(2) *(1 - z * x(2))) * n * (x(1) - n));
k3 = @(n,x) (x(3) * x(2) * (2 * n * x(1) - n)^2 * ( 1 - z * x(2)) + G * x(1) * (x(1) - 2 *n)* (x(1) - n) + x(4)^2 * G)./ (( G * x(3) - (x(1) - n)^2 * x(2) *(1 - z * x(2))) * n * (x(1) - n));
k4 = @(n,x) (x(4) * (x(1) + n)) ./ (n * (x(1)- n));
k5 = @(n,x) - x(5) ./ (n * (x(1)- n));

f = @(n,x) [k1(n,x);  k2(n,x);  k3(n,x);  k4(n,x); k5(n,x)];
[n,xa] = ode45(f,[0 1],[1-b 1/b 1-b 0.01 0.02]);

This runs for my install of Matlab.

However, the output is all NaNs since the function produce Infs; I may have introduced errors by replacing the brackets, but I don't know what the actual equations are so I'm going to leave that to you. :)

Answer 2

In general, you consult first the documentation. If not available, use the search engine of you choice, for instance https://www.google.de/search?q=matlab+ode45+example to find as first result https://de.mathworks.com/help/matlab/ref/ode45.html

There you find a multidimensional example

function dy = rigid(t,y)
    dy = zeros(3,1);    % a column vector
    dy(1) = y(2) * y(3);
    dy(2) = -y(1) * y(3);
    dy(3) = -0.51 * y(1) * y(2);
end

options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);
[T,Y] = ode45(@rigid,[0 12],[0 1 1],options);

plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')

that could be used as a blueprint. In the related page https://de.mathworks.com/help/matlab/math/ordinary-differential-equations.html they use a different syntax closer to yours for the van der Pol Equation in the first example

function dydt = vdp1(t,y)
    dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];
end

[t,y] = ode45(@vdp1,[0 20],[2; 0]);

plot(t,y(:,1),'-',t,y(:,2),'--')
title('Solution of van der Pol Equation, \mu = 1');
xlabel('time t');
ylabel('solution y');
legend('y_1','y_2')

Following these examples, rewrite your code as

G = 1.4;
g = 1.4;
k = 0;
z = 0;
b = 0.166667;



function dotx = dxdn(n,x) 
    t1 = n*(x1-n)    
    t2 = x(2)*(1 - z * x(2))

    den123 = ( G * x(3) - (x(1) - n)^2 * t2)

    k1 = ( x(1)*t1*t2 - 2 * x(3) * n *(1 - z * x(2)) - x(4)^2 * t2- G *x(3) *x(2) ) / (den123*n);   
    k2 = ( t2*(x(1) * x(2) * ( x(1) - 2 *n)*( x(1) - n) + 2* x(3) * n + x(4)^2 * x(2) ) ) / (den123 * t1);
    k3 = ( x(3) *  (2 * n * x(1) - n)^2 * t2 + G * x(1) * (x(1) - 2 *n)* (x(1) - n) + x(4)^2 * G ) / (den123*t1);
    k4 = ( x(4) * ( x(1) + n) ) / t1;
    k5 = - x(5) / t1;

    dotx = [k1;  k2;  k3;  k4; k5];
end 

[n,xa] = ode45(@dxdn,[0.001 1],[1-b; 1/b; 1-b; 0.01; 0.02]);

Because of the division by zero at n=0 you can not start the iteration in that singularity. In the code above, this is mitigated by starting with some (very) small positive n, you can also try starting with n=1e-8 or smaller. The slope will be very large in all components, so the integration may be slow, and the result might be not overly exact close to zero. For the correct handling of singular ODE ask in the math.stackexchange forum.