I have three coupled differential equation and need help to solve them

Initial condition: x(0) = 0.5, y(0)=0.9, z(0) = 7.5.

The problem is, I don't know how to call in Eq. 1. To solve these equations I am using ode45.

Kindly help me to solve this issue.

Thank you

[Star Strider]

Because of the presence of t in the denominator, the system is not defined (appoaches ∞) at .

The system likely has no finite solution.

Answer 1

[Sam Chak]

Yes @Captain Rituraj Singh,

The equations are rearranged. However, to simulate at , where division-by-zero occurs, I'm unable to comply.

syms x(t) y(t) z(t)
eqn1 = diff(x) == - (x/(4*t) + ((x^3)*z)/t + diff(z)) - y/t;
eqn2 = diff(y) == - (x*y/4 + (5*y/x)*diff(x)) - (x^4)/t;
eqn3 = diff(z) == - ((pi^2)/15)*(x/t + ((z^2)/x)*diff(x) - y/x);
eqns = [eqn1 eqn2 eqn3];
[V, S] = odeToVectorField(eqns) % convert to state equations for ode45 solver

V = 

S = 

Looking at the denominators, if , and , then singularities will occur at some points along this parabola of :

% Y1 = z(t), Y2 = x(t), Y3 = y(t)
sigma1 = @(x, z) 5926499560405365*z.^2 - 9007199254740992*x;
fimplicit(sigma1, [0 40 -8 8], 'LineWidth', 1.5)
grid on
xlabel({'$x(t)$'}, 'Interpreter', 'latex')
ylabel({'$z(t)$'}, 'Interpreter', 'latex')

Since the system cannot be simulated exactly at , and , then the initial sec is selected

odefcn = matlabFunction(V, 'vars', {'t', 'Y'})

odefcn = function_handle with value:
    @(t,Y)[((Y(2).^2.*-4.0+Y(1).^2.*Y(2)+Y(1).^2.*Y(3).*4.0+t.*Y(3).*4.0+Y(1).^3.*Y(2).^3.*4.0).*(-1.481624890101341e+15))./(t.*(Y(1).^2.*5.926499560405365e+15-Y(2).*9.007199254740992e+15));(Y(2).^2.*-3.674699746720117e+15+Y(1).*Y(2).^4.*9.007199254740992e+15+t.*Y(3).*5.926499560405365e+15+Y(2).*Y(3).*9.007199254740992e+15)./(t.*(Y(1).^2.*5.926499560405365e+15-Y(2).*9.007199254740992e+15));((Y(2).^6.*-3.602879701896397e+16-Y(2).^2.*Y(3).*7.349399493440234e+16+Y(2).*Y(3).^2.*1.801439850948198e+17+Y(1).^2.*Y(2).^5.*2.370599824162146e+16+t.*Y(3).^2.*1.185299912081073e+17-t.*Y(2).^3.*Y(3).*9.007199254740992e+15+Y(1).*Y(2).^4.*Y(3).*1.801439850948198e+17+t.*Y(1).^2.*Y(2).^2.*Y(3).*5.926499560405365e+15).*(-1.0./4.0))./(t.*(Y(1).^2.*5.926499560405365e+15-Y(2).*9.007199254740992e+15).*Y(2))]

tspan  = [1e-2 1.5];            % time span of simulation
y0     = [7.5 0.5 0.9];         % initial values: assumes z(1) = 7.5, x(1) = 0.5, y(1) = 0.9
opt    = odeset('RelTol', 1e-4, 'AbsTol', 1e-6);
[t, Y] = ode45(@(t, Y) odefcn(t, Y), tspan, y0);
plot(t, Y, 'LineWidth', 1.25), grid on, 
xlabel({'$t$'}, 'Interpreter', 'latex')

[Captain Rituraj Singh]

Hello Sam,

Thank you so much for the help. I tried a solving 3 coupled differential equation using your way but then it coming many warning and some errors as well: Please help me out to solve this issue. Here is the the my code

syms x(t) y(t) z(t) 


  a = 5.2037;
  b = 8.0260;
  c = 0.04829;
  d = 16.46;
 inv = 0.1973;

 tau0 = 0.1;%initial condition for t 
tauf = 10.0;

x0 = 0.50; %initial condition for x

s0 = c + (d*x0^3);

y0 = s0/(3*pi*tau0);% initial condition for y

z0 = 0.3; % initial condition for y
 %-------------------------------------------------------------
 %%%% Constant A, B and C
 %-------------------------------------------------------------
 

    q1 = 0.1973*3*(x^2)*z*cosh((0.1973*z)/(2*x));
    q2 = (0.1973*0.1973*(z^2)*x*sinh((0.1973*z)/(2*x)))/2;
    A = q1 - q2;
    
    q3 = (x.^3)*cosh((0.1973*z)/(2*x));
    q4 = (0.1973*(x)*z*sinh((0.1973*z)/(2*x)))/2;
    B = q3 + q4;
 
    C = (12*a*x.^(3)) + (3*A/(2*(pi^2)));
    
     w1 = (0.1973*2*(x^2)*z*cosh((0.1973*z)/(2*x)))/((pi^2)*(t));
     w2 = 4*a*(x.^4)/t;
     w3 = (3*B/(2*pi*pi))*0.1973*(diff(z));
     w4 = 0.1973*y/(t);

eqn1 = diff(x) == -(w1 + w2 + w3 - w4)/C;


    f1 = 2*a*x*y/(0.1973*3*b);
    f2 = y/(2*t);
    f3 = ((5*y)/(2*x))*diff(x);
    f4 = 8*a*(x.^4.0)/(0.1973*9*t);

eqn2 = diff(y) == -f1 - f2 + f3 + f4;

    r1 = -(pi^2)/(2*0.1973*B);
    r2 = (4*x/(3*t))*((c + d*(x.^3)) + (2*0.1973*(x^2)*z*cosh((0.1973*z)/(2*x))/(pi^2)) - ((0.1973*y)./x));
    r3 = (c + d*(x.^3) + 3*d*(x.^3) + (2*A/(pi^2)))*diff(x);    

eqn3 = diff(z) ==  r1*(r2 + r3);


eqns = [eqn1 eqn3 eqn3];

[V, S] = odeToVectorField(eqns);



odefcn = matlabFunction(V,'vars',{'t','V'});

Error using sym/matlabFunction>checkVarsSubset
Free variable 'Y' must be included in 'Vars' value.

Error in sym/matlabFunction>checkVars (line 240)
checkVarsSubset(vexpanded,funvars);

Error in sym/matlabFunction (line 158)
vars = checkVars(funvars,opts);

tspan = [tau0 tauf];

D     = [x0 y0 z0];

[t, Y] = ode45(@(t, Y) odefcn(t,V), tspan, D);

[Sam Chak]

Dear @Captain Rituraj Singh

I have made some corrections to a few lines. Consider voting 👍 my Answer above as a small token of appreciation.

syms x(t) y(t) z(t) 
a    = 5.2037;
b    = 8.0260;
c    = 0.04829;
d    = 16.46;
inv  = 0.1973;
tau0 = 0.1;             % initial condition for t 
tauf = 10.0;
x0   = 0.50;              % initial condition for x
s0   = c + (d*x0^3);
y0   = s0/(3*pi*tau0);    % initial condition for y
z0   = 0.3;               % initial condition for z

%-------------------------------------------------------------
% Constant A, B and C
%-------------------------------------------------------------
 
q1 = 0.1973*3*(x^2)*z*cosh((0.1973*z)/(2*x));
q2 = (0.1973*0.1973*(z^2)*x*sinh((0.1973*z)/(2*x)))/2;
A  = q1 - q2;
q3 = (x.^3)*cosh((0.1973*z)/(2*x));
q4 = (0.1973*(x)*z*sinh((0.1973*z)/(2*x)))/2;
B  = q3 + q4;
C  = (12*a*x.^(3)) + (3*A/(2*(pi^2)));

w1   = (0.1973*2*(x^2)*z*cosh((0.1973*z)/(2*x)))/((pi^2)*(t));
w2   = 4*a*(x.^4)/t;
w3   = (3*B/(2*pi*pi))*0.1973*(diff(z));
w4   = 0.1973*y/t;
eqn1 = diff(x) == -(w1 + w2 + w3 - w4)/C;

f1   = 2*a*x*y/(0.1973*3*b);
f2   = y/(2*t);
f3   = ((5*y)/(2*x))*diff(x);
f4   = 8*a*(x.^4.0)/(0.1973*9*t);
eqn2 = diff(y) == -f1 - f2 + f3 + f4;

r1   = -(pi^2)/(2*0.1973*B);
r2   = (4*x/(3*t))*((c + d*(x.^3)) + (2*0.1973*(x^2)*z*cosh((0.1973*z)/(2*x))/(pi^2)) - ((0.1973*y)./x));
r3   = (c + d*(x.^3) + 3*d*(x.^3) + (2*A/(pi^2)))*diff(x);    
eqn3 = diff(z) ==  r1*(r2 + r3);
eqns = [eqn1 eqn2 eqn3];                            % Correction in this line 1
[V, S] = odeToVectorField(eqns)

V = 

S = 

odefcn = matlabFunction(V, 'vars', {'t', 'Y'})      % Correction in this line #2

odefcn = function_handle with value:
    @(t,Y)[((Y(2).^4.*7.676280893357474e+51+Y(2).^7.*9.400276434808805e+37-Y(2).^3.*Y(3).*8.748940170530207e+51+Y(3).*2.425405993586982e+49-sinh((Y(1).*9.865e-2)./Y(2)).*Y(1).^2.*Y(2).^2.*4.848548970026535e+47+sinh((Y(1).*9.865e-2)./Y(2)).*Y(1).^2.*Y(2).^5.*4.372429292746232e+49-cosh((Y(1).*9.865e-2)./Y(2)).^2.*Y(1).^2.*Y(2).^4.*1.220776278130278e+49+cosh((Y(1).*9.865e-2)./Y(2)).^2.*Y(1).^2.*Y(2).^5.*1.220776278130278e+49-cosh((Y(1).*9.865e-2)./Y(2)).*Y(1).*Y(2).^2.*4.914900121669067e+48+cosh((Y(1).*9.865e-2)./Y(2)).*Y(1).*Y(2).^3.*1.47447003650072e+49-cosh((Y(1).*9.865e-2)./Y(2)).*Y(1).*Y(2).^5.*6.701118741161553e+51+cosh((Y(1).*9.865e-2)./Y(2)).*Y(1).*Y(2).^6.*7.144345228008612e+51+cosh((Y(1).*9.865e-2)./Y(2)).*sinh((Y(1).*9.865e-2)./Y(2)).*Y(1).^3.*Y(2).^3.*4.014319327918399e+47-cosh((Y(1).*9.865e-2)./Y(2)).*sinh((Y(1).*9.865e-2)./Y(2)).*Y(1).^3.*Y(2).^4.*4.014319327918399e+47).*(-2.160589846829774e-29))./(t.*(cosh((Y(1).*9.865e-2)./Y(2)).*Y(2).^2.*2.0e+4+sinh((Y(1).*9.865e-2)./Y(2)).*Y(1).*1.973e+3).*Y(2).*(Y(2).^3.*1.436445970991678e+18-3.982156237897728e+15));(Y(2).*(Y(2).^3.*2.419567526779995e+20+cosh((Y(1).*9.865e-2)./Y(2)).*Y(1).*Y(2).*2.221400516200497e+18-cosh((Y(1).*9.865e-2)./Y(2)).*Y(1).*Y(2).^2.*2.221400516200497e+18-2.683037472862266e+18).*(-1.97892936801441e-3))./(t.*(Y(2).^3.*1.436445970991678e+18-3.982156237897728e+15));((Y(2).^4.*2.05270293414876e+54-Y(2).^7.*7.404523286503226e+56+Y(2).^3.*Y(3).*4.211174100429751e+55-Y(3).*3.356371660383594e+53-t.*Y(2).*Y(3).*1.918174932234699e+53+t.*Y(2).^4.*Y(3).*6.919253008818116e+55+cosh((Y(1).*9.865e-2)./Y(2)).*Y(1).*Y(2).*Y(3).*2.41641951944048e+53-cosh((Y(1).*9.865e-2)./Y(2)).*Y(1).*Y(2).^2.*Y(3).*2.41641951944048e+53).*(-4.548045863171787e-38))./(t.*(Y(2).^3.*1.436445970991678e+18-3.982156237897728e+15))]

tspan  = [tau0 tauf];
D      = [z0 x0 y0];                                % Correction in this line #3
[t, Y] = ode45(@(t, Y) odefcn(t, Y), tspan, D);     % Correction in this line #4
plot(t, Y), grid on, xlabel('t')                    % Plot the result

[Captain Rituraj Singh]

Dear @Sam Chak,

Thank you so much for the helping me. Please, can you tell me how to extract the value of x, y, and z from Y and store these values in a file against t. I shall be very grateful to you.

Thank you once again.

Answer 2

[Sam Chak]

Hi @Captain Rituraj Singh

After this line

[t, Y] = ode45(@(t, Y) odefcn(t, Y), tspan, D);

the ode45 solver will generate the solutions in the Y array with respect ro the time vector t. To extract the solutions in the form that you are familiar with, you can do this:

z = Y(:, 1);    % because the odeToVectorField arranges them this way
x = Y(:, 2);
y = Y(:, 3);

and these solution vectors should appear on the Workspace. You can plot them if you wish

subplot(3,1,1)
plot(t, x)
subplot(3,1,2)
plot(t, y)
subplot(3,1,3)
plot(t, z)

If you want to keep and store the data {t, x, y, z} in a mat-file, where you can load and access it later, then delete the rest in the Workspace and enter this command:

save Captain.mat

If you want to load the data, enter this:

load Captain.mat

If you like the support, consider voting 👍 this Answer above as a small token of appreciation.

[Captain Rituraj Singh]

Thanks a lot @Sam Chak, you are a great help :)

Answer 3

[Sulaymon Eshkabilov]

Your given exercise is a type of implicitly defined ODEs and thus, you should employ ode15i - see this doc.

Answer 4

[Torsten]

You have a linear system of equations in dx/dt, dy/dt and dz/dt.

Solve it for these variables. Then you can write your system explicitly as

dx/dt = f1(t,x,y,z)
dy/dt = f2(t,x,y,z)
dz/dt = f3(t,x,y,z)

and use ode45 to solve (maybe by starting at t=1e-6 instead of t=0).