System of differential equations with many input arguments.

Hello,

I am trying to solve a system of five differential equations and get some plots back. The equations I have to solve are the following:

Where all correspond to derivatives of specific functions which in order to make the code for matlab I substituted them in the above functions. The code I am using for that is the following:

% l stands for \tau, u for \Theta, v for a, o for b (from utility) and j for ξ
syms S(t) P(t) c(t) h(t) z(t) A b d e g k l m n o r t u v w Y 
eqn1 = diff(S,t) == A*(b*S)^(1-g)*(z^g) - v*z - c - b*S*l;
eqn2 = diff(P,t) == u*z - d*P + j*c + (1-b)*S;
eqn3 = diff(h,t) ==  m*((b*S)-h);
eqn4 = diff(c,t) == (((b*(A*(1-g)*b^(1-g)*S^(-g)*z^g - l -r)) - ((A*g*(b*S)^(1-g)*z^(g-1) - v) * (1 - b + j*(d+r)) * (1/u))) * (u/(u - j*(v-g*A*(b*S)^(1-g)*z^(g-1)))))  +  ((1/(n*c^(-n)*(P^(-e)*b*S*h^(-o))^(1-n))) * ((exp(r*t)*w*m*b + b*(c*P^(-e)*h^(-o))^(1-n)*S^(-n)*b^(1-n)) + (j*e*P^(-e*(1-n)-1)*(c*b*S*h^(-o))^(1-n)) - ((1-n)*e*c^(-n)*P^(-e*(1-n)-1)*(b*S*h^(-o))^(1-n)*(u*z+(1-b)*S+j*c-d*P)) + (b*(1-n)*c^(-n)*(P^(-e)*h^(-o))^(1-n)*b^(1-n)*S^(-n)*(A*(b*S)^(1-g)*z^g - c - v*z - l*b*S))));
eqn5 = diff(z,t) == (((A*g*(b*S)^(1-g)*z^(g-1)-v) / (A*g*(g-1)*(b*S)^(1-g)*z^(g-2))) * ( d + ((e*P^(-e*(1-n)-1)*(c*b*S*h^(-o))^(1-n)*u) / (c^(-n)*(P*b*S*h^(-o))^(1-n))) - ((e*P^(-e*(1-n)-1)*(c*b*S*h^(-o))^(1-n)*j) / (c^(-n)*(P^(-e)*b*S*h^(-o))^(1-n))) + b*(A*(1-g)*b^(1-g)*S^(-g)*z^g - l + ((w*m*(u-j*(v-A*g*(b*S)^(1-g)*z^(g-1)))) / (exp(-r*t)*u*c^(-n)*(P^(-e)*b*S*h^(-o))^(1-n)))) + ((b^(1-n)*S^(-n)*(c*P^(-e)*h^(-o))^(1-n) * (u-j*(v-A*g*(b*S)^(1-g)*z^(g-1)))) / (u*e^(-n)*(P^(-e)*b*S*h^(-o))^(1-n))) + ((1-b)*(v-A*g*(b*S)^(1-g)*z^(g-1))*(1/u)))) - ((b*z/S) * (A*(b*S)^(1-g)*z^g - v*z - c - b*S*l));

[VF,Subs] = odeToVectorField(eqn1,eqn2,eqn3,eqn4,eqn5);
odefcn = matlabFunction(VF, 'Vars',{t,A,Y,b,d,e,g,k,l,m,n,o,r,u,v,w});
A=3;
b=0.6;
d=0.3;
e=0.4;
g=0.3;
k=4;
l=2;
m=0.7;
n=2;
o=0.6;
r=0.5;
u=10;	
v=2;
w=0.7;

tspan = [0 50];
y0 = [100 60 40 30 50];
[t,Y] = ode45(@(t,Y) odefcn(t,A,Y,b,d,e,g,k,l,m,n,o,r,t,u,v,w) , tspan, y0);
P = Y(:,1);	
S = Y(:,2);
h = Y(:,3);
c = Y(:,4);
z = Y(:,5);
x = b*S;		
q = A*((b*S).^(1-g)).*(z.^g);

figure(1)
plot(c,P)
title(sprintf('c-P'))
xlabel('Consumption (c)')
ylabel('Pollution (P)')

figure(2)
plot(c,x)
title(sprintf('C-x'))
xlabel('Consumption (c)')
ylabel('Recycling (x)')

figure(3)
plot(q,P)
title(sprintf('q-P'))
xlabel('Output (q)')
ylabel('Pollution (P)')

figure(4)
plot(q,x)
title(sprintf('q-x'))
xlabel('Output (q)')
ylabel('Recycling (x)')

I've used that code before for different functions and it was working perfectly but with those differential equations it doesn't work and it give me back the following error:

Error using
symengine>@(t,A,Y,b,d,e,g,k,l,m,n,o,r,u,v,w)[-d.*Y(1)+u.*Y(5)+Y(4).*1i-Y(2).*(b-1.0);-v.*Y(5)-Y(4)+A.*(b.*Y(2)).^(-g+1.0).*Y(5).^g-b.*l.*Y(2);m.*(b.*Y(2)-Y(3));-(u.*(b.*(l+r+A.*b.^(-g+1.0).*Y(2).^(-g).*Y(5).^g.*(g-1.0))-((v-A.*g.*(b.*Y(2)).^(-g+1.0).*Y(5).^(g-1.0)).*(-b+d.*1i+r.*1i+1.0))./u))./(u-v.*1i+A.*g.*(b.*Y(2)).^(-g+1.0).*Y(5).^(g-1.0).*1i)+(Y(4).^n.*(b.*Y(1).^(-e).*Y(2).*Y(3).^(-o)).^(n-1.0).*(e.*Y(1).^(e.*(n-1.0)-1.0).*(b.*Y(2).*Y(3).^(-o).*Y(4)).^(-n+1.0).*1i+b.*m.*w.*exp(r.*t)+b.*b.^(-n+1.0).*Y(2).^(-n).*(Y(1).^(-e).*Y(3).^(-o).*Y(4)).^(-n+1.0)-e.*Y(1).^(e.*(n-1.0)-1.0).*Y(4).^(-n).*(n-1.0).*(b.*Y(2).*Y(3).^(-o)).^(-n+1.0).*(d.*Y(1)-u.*Y(5)-Y(4).*1i+Y(2).*(b-1.0))+b.*b.^(-n+1.0).*(Y(1).^(-e).*Y(3).^(-o)).^(-n+1.0).*Y(2).^(-n).*Y(4).^(-n).*(n-1.0).*(v.*Y(5)+Y(4)-A.*(b.*Y(2)).^(-g+1.0).*Y(5).^g+b.*l.*Y(2))))./n;(b.*Y(5).*(v.*Y(5)+Y(4)-A.*(b.*Y(2)).^(-g+1.0).*Y(5).^g+b.*l.*Y(2)))./Y(2)-((b.*Y(2)).^(g-1.0).*Y(5).^(-g+2.0).*(v-A.*g.*(b.*Y(2)).^(-g+1.0).*Y(5).^(g-1.0)).*(d-b.*(l+A.*b.^(-g+1.0).*Y(2).^(-g).*Y(5).^g.*(g-1.0)-(m.*w.*exp(r.*t).*Y(4).^n.*(u-v.*1i+A.*g.*(b.*Y(2)).^(-g+1.0).*Y(5).^(g-1.0).*1i).*(b.*Y(1).^(-e).*Y(2).*Y(3).^(-o)).^(n-1.0))./u)-((v-A.*g.*(b.*Y(2)).^(-g+1.0).*Y(5).^(g-1.0)).*(b-1.0))./u-e.*Y(1).^(e.*(n-1.0)-1.0).*Y(4).^n.*(b.*Y(2).*Y(3).^(-o).*Y(4)).^(-n+1.0).*(b.*Y(1).^(-e).*Y(2).*Y(3).^(-o)).^(n-1.0).*1i+e.*u.*Y(1).^(e.*(n-1.0)-1.0).*Y(4).^n.*(b.*Y(1).*Y(2).*Y(3).^(-o)).^(n-1.0).*(b.*Y(2).*Y(3).^(-o).*Y(4)).^(-n+1.0)+(b.^(-n+1.0).*e.^n.*Y(2).^(-n).*(u-v.*1i+A.*g.*(b.*Y(2)).^(-g+1.0).*Y(5).^(g-1.0).*1i).*(Y(1).^(-e).*Y(3).^(-o).*Y(4)).^(-n+1.0).*(b.*Y(1).^(-e).*Y(2).*Y(3).^(-o)).^(n-1.0))./u))./(A.*g.*(g-1.0))]
Too many input arguments.

Error in @(t,Y)odefcn(t,A,Y,b,d,e,g,k,l,m,n,o,r,t,u,v,w)

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

So now I do not know how to move on. I understand that the system of differential equations I am using has many arguments as the error message says as well, but most of them are parameters which I define giving numerical values for them. Shouldn't Matlab be able to solve that?

Any help would be very much appreciated (any idea, thought, solution or anything else).

Thank you so much in advance.

Answer 1

[Stephan]

Try:

% l stands for \tau, u for \Theta, v for a, o for b (from utility) and j for ξ
syms S(t) P(t) c(t) h(t) z(t) A b d e g k l m n o r t u v w Y 
eqn1 = diff(S,t) == A*(b*S)^(1-g)*(z^g) - v*z - c - b*S*l;
eqn2 = diff(P,t) == u*z - d*P + j*c + (1-b)*S;
eqn3 = diff(h,t) ==  m*((b*S)-h);
eqn4 = diff(c,t) == (((b*(A*(1-g)*b^(1-g)*S^(-g)*z^g - l -r)) - ((A*g*(b*S)^(1-g)*z^(g-1) - v) * (1 - b + j*(d+r)) * (1/u))) * (u/(u - j*(v-g*A*(b*S)^(1-g)*z^(g-1)))))  +  ((1/(n*c^(-n)*(P^(-e)*b*S*h^(-o))^(1-n))) * ((exp(r*t)*w*m*b + b*(c*P^(-e)*h^(-o))^(1-n)*S^(-n)*b^(1-n)) + (j*e*P^(-e*(1-n)-1)*(c*b*S*h^(-o))^(1-n)) - ((1-n)*e*c^(-n)*P^(-e*(1-n)-1)*(b*S*h^(-o))^(1-n)*(u*z+(1-b)*S+j*c-d*P)) + (b*(1-n)*c^(-n)*(P^(-e)*h^(-o))^(1-n)*b^(1-n)*S^(-n)*(A*(b*S)^(1-g)*z^g - c - v*z - l*b*S))));
eqn5 = diff(z,t) == (((A*g*(b*S)^(1-g)*z^(g-1)-v) / (A*g*(g-1)*(b*S)^(1-g)*z^(g-2))) * ( d + ((e*P^(-e*(1-n)-1)*(c*b*S*h^(-o))^(1-n)*u) / (c^(-n)*(P*b*S*h^(-o))^(1-n))) - ((e*P^(-e*(1-n)-1)*(c*b*S*h^(-o))^(1-n)*j) / (c^(-n)*(P^(-e)*b*S*h^(-o))^(1-n))) + b*(A*(1-g)*b^(1-g)*S^(-g)*z^g - l + ((w*m*(u-j*(v-A*g*(b*S)^(1-g)*z^(g-1)))) / (exp(-r*t)*u*c^(-n)*(P^(-e)*b*S*h^(-o))^(1-n)))) + ((b^(1-n)*S^(-n)*(c*P^(-e)*h^(-o))^(1-n) * (u-j*(v-A*g*(b*S)^(1-g)*z^(g-1)))) / (u*e^(-n)*(P^(-e)*b*S*h^(-o))^(1-n))) + ((1-b)*(v-A*g*(b*S)^(1-g)*z^(g-1))*(1/u)))) - ((b*z/S) * (A*(b*S)^(1-g)*z^g - v*z - c - b*S*l));

[VF,Subs] = odeToVectorField(eqn1,eqn2,eqn3,eqn4,eqn5);
odefcn = matlabFunction(VF, 'Vars',{t,Y,A,b,d,e,g,k,l,m,n,o,r,u,v,w});
A=3;
b=0.6;
d=0.3;
e=0.4;
g=0.3;
k=4;
l=2;
m=0.7;
n=2;
o=0.6;
r=0.5;
u=10;	
v=2;
w=0.7;

tspan = [0 50];
y0 = [100 60 40 30 50];
[t,Y] = ode45(@(t,Y) odefcn(t,Y,A,b,d,e,g,k,l,m,n,o,r,u,v,w) , tspan, y0);
P = Y(:,1);	
S = Y(:,2);
h = Y(:,3);
c = Y(:,4);
z = Y(:,5);
x = b*S;		
q = A*((b*S).^(1-g)).*(z.^g);

figure(1)
plot(c,P)
title(sprintf('c-P'))
xlabel('Consumption (c)')
ylabel('Pollution (P)')

figure(2)
plot(c,x)
title(sprintf('C-x'))
xlabel('Consumption (c)')
ylabel('Recycling (x)')

figure(3)
plot(q,P)
title(sprintf('q-P'))
xlabel('Output (q)')
ylabel('Pollution (P)')

figure(4)
plot(q,x)
title(sprintf('q-x'))
xlabel('Output (q)')
ylabel('Recycling (x)')

3 reasons:

[t,Y] = ode45(@(t,Y) odefcn(t,A,Y,b,d,e,g,k,l,m,n,o,r,t,u,v,w) , tspan, y0);
%                             ^ ^                     ^
%                             | |                     |
%                             -----> Change order     ----> t already at first position

odefcn = matlabFunction(VF, 'Vars',{t,A,Y,b,d,e,g,k,l,m,n,o,r,u,v,w});
%                                   ^ ^                     
%                                   | |                     
%                  SAME HERE        -----> Change order    

[Myrto Kasioumi]

Thak you so much. Now it doesn't give me any error back but still it doesn't give me results (it runs for many minutes and still nothing). Do you think it's because of the values I am using or is it because of the form of the functions? Do you have any suggestions for that problem too?

Thanks in advance!

[Stephan]

Your system appears to be stiff. try to use ode15s or similar stiff solvers. Also you could try to use odeset and try to reduce tolerance values if this is possible for your case. Do you need to calculate the full time between 0 and 50? This could also help. It is a bit of experimenting.

[Myrto Kasioumi]

Thank you that helped me a lot.