Explicit solution not found

I have a large system of equations. I want to look at their derivatives.

syms X y1(X) y2(X) y3(X) y4(X) y5(X) y6(X) y7(X) y8(X) y9(X) y10(X) y11(X) y12(X) y13(X) y14(X) y15(X) y16(X) y17(X) y18(X) y19(X) y20(X) y21(X) y22(X) y23(X)

k1i = 0.001; k1 = 0.01; k2i = 0.05; k2ii = 1; k2iii = 1;
k3i=0;k3ii=10;k4i=40;k4=40;J3=0.01;J4=0.01;
k7=1;k8=0.2;J7=0.01;J8=0.01;
k9=0.1;k10=0.04;J9=0.01;J10=0.01;
k11=0.2;k12=0.2;k12i=2;k12ii=100;
lp=1000;lm=1;kssim=100;kdsimi=0.0001;kdsim=10;
k13i=0;k13ii=0.05;k14i=0.02;k14=1;
k15i=0.25;k15ii=0;k16i=0.01;k16ii=2;J15=0.1;J16=0.1;
k21i=0;k21=10;k22=0.5;Jtfb=1;
lcp=40;lcm=1;mm=0.1;mp=0.1;
k1111=0.1;js=0.1;js2=0.1;ks=0.1;
ks1=0.1;ks2=0.1;sm=0.1;sp=0.1; kscdh=0.1; vdcdh=0.1; kstfb=0.1;vdtfb=0.1;

eq1=diff(y1)==k1i+k1.*y9-y13.*y1+mp.*y19; 
eq2=diff(y2)==y19-y3;
eq3=diff(y3)==mm.*y18.*y1 -mp.*y3 -lp.*y3.*(y4-y2-y23)+lm.*(y19-y3);
eq4=diff(y4)==(ks+ks1.*y12).*(y20-y4)./(js+y20-y4)-(ks2.*y4)./(js2+y4);
eq5=diff(y5)==(k3i+k3ii.*y7).*(y21-y5)./(J3+y21-y5)-(k4i.*y12+k4.*y3).*y5./(J4+y5);
eq6=diff(y6)==k9.*(1-y6).*y3./(J9+1-y6)-k10.*y6./(J10+y6);
eq7=diff(y7)==k7.*y6.*(1-y7)./(J7+1-y7)-k8.*y7./(J8+y7);
eq8=diff(y8)==k11-y14.*y8;
eq9=diff(y9)==(k21i+k21.*y3).*(y22-y9)./(Jtfb+y22-y9)-k22.*y9./(Jtfb+y9);
eq10=diff(y10)==k15i.*(1-y10)./(J15+1-y10)-(k16i+k16ii.*y12).*y10./(J16+y10);
eq11=diff(y11)==k13i+k13ii.*y10-(k14i+k14.*y3).*y11;
eq12=diff(y12)==k13i+k13ii.*y10+lcm.*y16-lcp.*y12.*y17-(k14i+k14.*y3).*y12+sm.*y23-sp.*y12.*(y4-y2-y23);
eq13=diff(y13)==k2i+k2ii.*y5+k2iii.*y7;
eq14=diff(y14)==k12+k12i.*y12+k12ii.*y3;
eq15=diff(y15)==kdsimi+kdsim.*y3;
eq16=diff(y16)==y11-y12-y23;
eq17=diff(y17)==y8-y16;
eq18=diff(y18)==k1i+k1.*y9-k1111.*y18+mp.*y19;
eq19=diff(y19)== mm.*y18.*y1-mp.*y19;
eq20=diff(y20)==kssim-y15.*y20;
eq21=diff(y21)==kscdh-vdcdh.*y21;
eq22=diff(y22)==kssim-y15.*y22;
eq23=diff(y23)==-sm.*y23+sp.*y12.*(y4-y2-y23);

q=dsolve(eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10,eq11,eq12,eq13,eq14,eq15,eq16,eq17,eq18,eq19,eq20,eq21,eq22,eq23,'X');

Can anyone tell me why it can't find a solution? I had a smaller version of this, and it also wouldn't evaluate. I would like to evaluate the derivatives at one or two points to find the slopes. I'm not 100% on how the symbolic toolbox works, so can someone shed some light on this? Any info would be greatly appreciated.

Answer 1

[Walter Roberson]

syms x y1 y2 y3 y4 y5 y6 y7 y8 y9 y10 y11 y12 y13 y14 y15 y16 y17 y18 y19 y20 y21 y22 y23

and

syms Jy9 kdy4i kdy4 ksy4 ksy9 vdy9

The second list is needed to account for undefined variables in your code. Perhaps some of those are typing mistakes in your code. For example, perhaps kdy4i is intended to be k*diff(y4(X),X)*1i -- we don't know.

You should be replacing each of those diff() to diff() with respect to a variable, such as diff(y4(X),X)

[Renee]

Ah crap sorry. I switched the variable names with y() to try to simplify a bit, so that's where they are coming from. The code I'm actually running has the required definitions.

[Renee]

I also noticed I declared my

syms

differently. Fixed.

[Walter Roberson]

eq22 uses the undefined vdsim

[Renee]

So sorry, fixed.

[Walter Roberson]

Okay, that code could plausibly work. It is, though, a fair sized system, and it might take a lot of time and memory to resolve.

[Renee]

Do you think that's why its not finding a solution? dsolve itself only takes 2 seconds to return "solution not found". Are you supposed to write the solve function differently when they all have to be solved simultaneously?

[Walter Roberson]

I did not try the Symbolic Toolbox dsolve() on that system. I tried Maple's dsolve() on the system, and it quickly chewed through enough memory that I needed to kill the computation. I was expecting something like that might happen. Unfortunately it does not take much for a system of simultaneous differential equations to "explode".

My rule of thumb for nonlinear simultaneous non-differential equations with powers as low as 2, is that 7 simultaneous equations is about the break-point, where 6 low-degree non-differential equations can probably be solved in a reasonable time, and 8 probably cannot be solved in a reasonable time and memory, and 7 might or might not work. As the degree increases, the number of practical equations goes down. That said, 4 such equations is often the limit of "fast and compact", and when it starts getting to 5 equations or more it is not uncommon for me to need to start helping it, especially by choosing particular branches to focus on when multiple roots start popping up.

For simultaneous differential equations, I find that it is not at all uncommon for 2 simultaneous differential equations to either fail or be at the limit of practicality. The only time I have seen more work is if they were especially simple equations. I would be amazed to see 22 simultaneous differential equations resolved symbolically in a form that was comprehendable.

For that many equations, you are very likely need to switch to numeric solutions, which will require that you create numeric boundary conditions.

It is the non-numeric boundary conditions that seems to especially bog down symbolic resolution of simultaneous differential equations. "This boundary condition has to be the x value such that integral( this_complicated_function_with_some_unbound_variables(x), x in some_range) exactly equals one of the 207th roots of this 500 page expression".

[Renee]

Thanks for the info! I added in some initial conditions, "y(i)==0" for all variables and same result. AFAIK that's the extent of boundary conditions I can change, am I correct? So do you think its impossible then?

[Walter Roberson]

You can also control the boundary conditions for the derivatives.

[Walter Roberson]

The equations such as for y4, that involve division by the same function, are not easy to automatically deal with. If you try to solve one of them in isolation without any boundary conditions then you get an expression involving roots and arctan and ln, but that expression fails at X = 0 because it encounters ln(0). When you add in the boundary conditions then it looks to me as if the solutions might have to be that the function is a constant 0, but I am not certain of that.

[Renee]

Thanks for pointing that out, I didn't notice that. I changed the initial bounds to be nonzero, and I know that some equations' derivatives are 0 so I added the bounds "diff(yi(X))==0" for about half the equations. Now I am getting the error "Cannot reduce to the square system because the number of equations differs from the number of indeterminates." Am I too hopeful in thinking this is progress? I can try to find an additional equation to solve this perhaps.