System of differential equations, unable to find explicit solution with dsolve

Hi there!

I am trying to solve the following system of differential equations (4 equations and 4 unknowns):

where and and

So, when I substitute all the derivatives I get the following system of differential equations:

I don't have any initial conditions or any values for the parameters.

I am not familiar with Matlab but after searching online I wrote this code which doesn't work:

syms S(t) P(t) c(t) z(t) A b g v k u d r w	% k stands for τ, u for θ and v for a
eqn1 = diff(S,t) == A*(b*S)^(1-g)*(z^g) - k*z - c - b*S*(k+1)
eqn2 = diff(P,t) == u*z - d*P + (1-b)*S
eqn3 = diff(c,t) == (((1-b)/u) * (v - A*g*(b*S)^(1-g)*z^(g-1)) + b * (A*(1-g)*b^(-g)*S^(-g)*z^(g) -k -1) - r) * c
eqn4 = diff(z,t) == ((v - A*g*(b*S)^(1-g)*z^(g-1))/(A*(g-1)*g*(b*S)^(1-g)*z^(g-2))) * ((d - ((1-b) * (A*g*(b*S)^(1-g)*z^(g-1) - v))/u)) + (b / (A*(g-1)*g*(b*S)^(1-g)*z^(g-2))) * (((A*g*(b*S)^(1-g)*z^(g-1) - v)) * (A*b^(2-g)*(1-g)*S^(-g)*z^g - 1 - k) - (A*g*(1-g)*b^(1-g)*S^(-g)*z^(g-1)) * (A*(b*S)^(1-g)*z^g) - v*z - b*S*(k+1) - c) - (((1-w)*u *(1/P)) / (A*(g-1)*g*(b*S)^(1-g)*z^(g-2)))
sol = dsolve(eqn1,eqn2,eqn3,eqn4)

What I am getting back is the following error and an empty matrix of solutions:

Warning: Unable to find explicit solution.
> In dsolve (line 190) 

But I am sure that there is a solution for that system.

The version of Matlab I am using is: 9.7.0.1319299 (R2019b) Update 5.

Can't seem to fix this problem even though I tried many things and I don't understand what I am doing wrong. I would really appreciate any help.

Thanks in advance!

Answer 1

[Star Strider]

Not all differential equations (or integrals) have analytic sollutions. Yours are among them.

Try this:

% k stands for \tau, u for \Theta and v for a
syms S(t) P(t) c(t) z(t) A b g v k u d r w t Y 
eqn1 = diff(S,t) == A*(b*S)^(1-g)*(z^g) - k*z - c - b*S*(k+1);
eqn2 = diff(P,t) == u*z - d*P + (1-b)*S;
eqn3 = diff(c,t) == (((1-b)/u) * (v - A*g*(b*S)^(1-g)*z^(g-1)) + b * (A*(1-g)*b^(-g)*S^(-g)*z^(g) -k -1) - r) * c;
eqn4 = diff(z,t) == ((v - A*g*(b*S)^(1-g)*z^(g-1))/(A*(g-1)*g*(b*S)^(1-g)*z^(g-2))) * ((d - ((1-b) * (A*g*(b*S)^(1-g)*z^(g-1) - v))/u)) + (b / (A*(g-1)*g*(b*S)^(1-g)*z^(g-2))) * (((A*g*(b*S)^(1-g)*z^(g-1) - v)) * (A*b^(2-g)*(1-g)*S^(-g)*z^g - 1 - k) - (A*g*(1-g)*b^(1-g)*S^(-g)*z^(g-1)) * (A*(b*S)^(1-g)*z^g) - v*z - b*S*(k+1) - c) - (((1-w)*u *(1/P)) / (A*(g-1)*g*(b*S)^(1-g)*z^(g-2)));
[VF,Subs] = odeToVectorField(eqn1,eqn2,eqn3,eqn4)
odefcn = matlabFunction(VF, 'Vars',{t,A,Y,b,d,g,k,r,u,v,w})

The ‘Subs’ output is important because it tells how the ‘Y’ vector elements are assigned, so ‘Y(1)=Subs(1)’ and so for the rest.

Supply scalar values for the other variables, and (probably) use:

@(t,Y) odefcn(t,A,Y,b,d,g,k,r,u,v,w) 

as the function argument to the ODE solver of your choice.

.

[Myrto Kasioumi]

Thank you so much. That really helped me!

[Star Strider]

As always, my pleasure!

[Myrto Kasioumi]

I set values for all the parameters as follows:

A=3;
b=0.4;
g=0.3;
v=2;
k=4;
u=10;
d=0.2;
r=0.5;
w=0.7;

And I also set initial values for all the dependent variables that I have and time interval so the calculation of the ode function like so:

y0 = [100 50 60 30];
tspan = [0 50];

So the total code I have so far is:

% k stands for \tau, u for \Theta and v for a
syms S(t) P(t) c(t) z(t) A b g v k u d r w t Y 
eqn1 = diff(S,t) == A*(b*S)^(1-g)*(z^g) - k*z - c - b*S*(k+1);
eqn2 = diff(P,t) == u*z - d*P + (1-b)*S;
eqn3 = diff(c,t) == (((1-b)/u) * (v - A*g*(b*S)^(1-g)*z^(g-1)) + b * (A*(1-g)*b^(-g)*S^(-g)*z^(g) -k -1) - r) * c;
eqn4 = diff(z,t) == ((v - A*g*(b*S)^(1-g)*z^(g-1))/(A*(g-1)*g*(b*S)^(1-g)*z^(g-2))) * ((d - ((1-b) * (A*g*(b*S)^(1-g)*z^(g-1) - v))/u)) + (b / (A*(g-1)*g*(b*S)^(1-g)*z^(g-2))) * (((A*g*(b*S)^(1-g)*z^(g-1) - v)) * (A*b^(2-g)*(1-g)*S^(-g)*z^g - 1 - k) - (A*g*(1-g)*b^(1-g)*S^(-g)*z^(g-1)) * (A*(b*S)^(1-g)*z^g) - v*z - b*S*(k+1) - c) - (((1-w)*u *(1/P)) / (A*(g-1)*g*(b*S)^(1-g)*z^(g-2)));
[VF,Subs] = odeToVectorField(eqn1,eqn2,eqn3,eqn4);
odefcn = matlabFunction(VF, 'Vars',{t,A,Y,b,d,g,k,r,u,v,w});

A=3;
b=0.4;
g=0.3;
v=2;
k=4;
u=10;
d=0.2;
r=0.5;
w=0.7;

tspan = [0 50];
y0 = [100 50 60 30];
[t,Y] = ode45(@(t,Y) odefcn(t,A,Y,b,d,g,k,r,u,v,w) , tspan, y0)

The ode function works but I still get a warning back saying:

Warning: Failure at t=4.800543e-03.  Unable to meet integration tolerances without reducing the step
size below the smallest value allowed (1.387779e-17) at time t. 
> In ode45 (line 360)

I also used the dbstop if naninf command like that:

y0 = [100 50 60 30];
dbstop if naninf
[t,Y] = ode45(@(t,Y) odefcn(t,A,Y,b,d,g,k,r,u,v,w) , tspan, y0)

and what it gives me back is the following:

NaN/Inf breakpoint hit for ode45.m on line 256.
256     rh = norm(f0 ./ max(abs(y),threshold),inf) / (0.8 * rtol^pow);

Even though the code is running for some t, I don't get back values for the dependent variable Y which contains (P,S,c,z). What I get back is the values for the independent variable t and then for the dependent variables I get a matrix with four columns but the first 3 are all 0 and I get values only for the last one (I also changed the initial values of y0 and I managed to get values for the first column as well but it was always constant and still I didn't get values for the second and third column).

Then I tried to seperate the columns into different vectors so I can see the values clearer and make some plots like that:

Y1 = Y(1:end,1)
Y2 = Y(1:end,2)
Y3 = Y(1:end,3)
Y4 = Y(1:end,4)

but when I am doing that all the values change and I get different results than what I see in the initial matrix that ode gives me back.

So my problems are:

• the calculations stop at some t before maxt=50
• I don't get back values for all the dependent variables
• when I am trying to seperate the columns of the matrix of ode solution, the values change

Any help would be appreciated.

Thanks a lot.

[Star Strider]

The problem that throws the Warning:

Warning: Failure at t=4.800543e-03.  Unable to meet integration tolerances without reducing the step
size below the smallest value allowed (1.387779e-17) at time t. 

is that ode45 has encountered a singularity and cannot continue the integration. The only solution for that is to examine ‘odefcn’ carefully to see where the ±Inf values occur.

I assume those equations are from a published paper, so if the paper does not appear to have had problems integrating with those equations and initial conditions, it would be worthwhile to examine your symbolic code carefully to be sure the equations are correct. If they are, then see what solver the paper used. Some solvers (especially fixed-step solvers) are not as sophisticated as the MATLAB solvers, and so could integrate the system (possibly erroneously) to produce the published data. And there is always the possibility that the paper’s results are not reproducible. (You would not be the first to encounter that.)

I also got the Warning with different MATLAB ODE solvers. Looking at the output, ‘z’ seems to be the problem, so concentrate your efforts on it.

I used this to plot them:

sbp = size(Y,2);
figure
for k = 1:sbp
    subplot(sbp,1,k)
    if k < sbp
        plot(t,Y(:,k))
    else
        semilogy(t,Y(:,k))
    end
    grid
    title(string(Subs(k)))
end

.

[Myrto Kasioumi]

Thanks a lot for your help again!!

[Star Strider]

As always, my pleasure!

Answer 2

[Myrto Kasioumi]

Hello again!

I am still trying to solve the above problem but now I made it a little bit harder. So the basic problem is still the same.

I am trying to solve the following system of differential equations:

So I am using the following code:

% k stands for \tau, u for \Theta, v for a, h for \hta and A for Ao
syms S(t) P(t) c(t) z(t) A b g h v k l u d r w t Y 

eqn1 = diff(S,t) == (A*exp(h*t)*b*S)^(1-g) * z^g - v*z - c - A*exp(h*t)*b*S*(k*exp(h*t*(-1))+1);
eqn2 = diff(P,t) == u*z - d*P + (1-b)*S;
eqn3 = diff(c,t) == ( ((A*exp(h*t)*b) * ((1-g)*b*z^g*(A*exp(h*t)*S*b)^(-g)*(A*exp(h*t)+S) - 1 - k*exp(l*t*(-1))))  + (((1-b)/u) * (v - g*z^(g-1)*(A*exp(h*t)*S*b)^(1-g))) - r ) * c ;
eqn4 = diff(z,t) == ( ( (g*(A*exp(h*t)*S*b)^(1-g)*z^(g-1) - v)/(g*(g-1)*(A*exp(h*t)*S*b)^(1-g)*z^(g-2)) ) * ( d + ((c*u*(1-w))/(P*w*(v-g*(A*exp(h*t)*S*b)^(1-g)*z^(g-1)))) + ((A*exp(h*t)*b) * ((1-g)*z^g*(A*exp(h*t)*b*S)^(-g)*(A*exp(h*t)*b + b*S) - 1 - (k*exp(l*t*(-1))))) + (((1-b)*(v-g*(A*exp(h*t)*b*S)^(1-g)*z^(g-1)))/(u)) ) ) + ( ( (b*(A*exp(h*t)+S)*z)/(A*exp(h*t)*S) ) * ( S*A*exp(h*t)*h + (A*exp(h*t)) * ((A*exp(h*t)*b*S)^(1-g)*z^g - v*z - c - (1+(k*exp(l*t*(-1))))*A*exp(h*t)*b*S) ) );

[VF,Subs] = odeToVectorField(eqn1,eqn2,eqn3,eqn4);
odefcn = matlabFunction(VF, 'Vars',{t,A,Y,b,d,g,h,k,l,r,u,v,w});

A=3;
b=0.4;
g=0.3;
h=0.3;
v=2;
k=4;
l=0.7;
u=10;
d=0.2;
r=0.5;
w=0.7;

tspan = [0 1000];
y0 = [100 60 30 50];    
[t,Y] = ode45(@(t,Y) odefcn(t,A,Y,b,d,g,h,k,l,r,u,v,w) , tspan, y0);

P_1 = Y(:,1);	
S_1 = Y(:,2);
c_1 = Y(:,3);
z_1 = Y(:,4);
AF_1 = A*exp(h*t);
x_1 = b*S.*AF;

Now I wanna see where the marginal product of x is lower than the MP of z so I am doing that:

MPx = (1-g)*(b*A*exp(h*t).*S_1).^(-g).*z_1.^g*b.*(A*exp(h*t)+S_1);
MPz = g*(b*A*exp(h*t).*S_1).^(1-g).*z_1.^(g-1);
m = find(MPx(:,1)<MPz(:,1));
i = size(m);

For the values that MPx is lower than MPz I will run the previous code again but for g = 0.7 now instead of 0.3 and then I wanna replace the values of all P_1, S_1, c_1, z_1, x_1 with the new values that I will get from the new code but only in the specific rows for which hold that MPx < MPz.

So I am running this part again:

g=0.7;
tspan = [0 1000];
y0 = [100 60 30 50];    %even with initial values equal to 1 for all, the plots are the same
[t,N] = ode45(@(t,N) odefcn(t,A,N,b,d,g,h,k,l,r,u,v,w) , tspan, y0);

P_2 = N(:,1);	%check the size by size(P)
S_2 = N(:,2);
c_2 = N(:,3);
z_2 = N(:,4);
AF_2 = A*exp(h*t);
x_2 = b*S.*AF;

And now I wanna do something that will help me substitutes the values of P_1 etc with the values of P_2 etc in the rows for which MPx < MPz.

Any help would be really appreciated.

Thanks in advance.

[Star Strider]

I am finding this difficult to follow.

It will likely be easier to do that substitution if tspan has the same defined length for both integrations, for example:

tspan = linspace(0, 1000, 250);

creating outputs with row lengths of 250 here (choose whatever length you want). Then the substitution, using either ‘m’ or a logical indexing vector for the rows, should be straightforward.

[Myrto Kasioumi]

Thanks a lot.

[Star Strider]

As always, my pleasure!