how to obtain matrix with two symbolic equations with dsolve

clc, clear
a1=28*7*400;                        %%ft^2
po=40.6;                            %%psi
dpio=9;                           %%psi
qf=1525;                            %%gpm
qo=qf;
pp=16.4;                            %%psi
qp=1234;                            %%gpm
depr1=24.9;                         %%psi
delp=depr1-pp;                      %%psi
delpr2=18.2;                        %%psi
PH=360;                                     %%ft
lp=2.76*.264/(60)^2*14.5;                   %%gal/m*min^2/psi
k2=1.75*10^-4*60^1.67*14.5/264.172^1.67;     %%psi/(gal/min)^1.67
k3=2.27*10^-2;                                  %%(m^3/h)^.6/m^2
g=32.17;                            %%ft/s^2
rhoh20=62.4;                        %%lbs/ft^3
pf=rhoh20*g*PH*1.488*.000145;       %%psi
dpo=pf-po;                     %%psi
eqns=ro_cp_model(lp,qo,dpio,k3,k2,dpo);

function dpq=ro_cp_model(lp,qo,dpio,k3,k2,dpo)

syms jjw dp(x) Q(x)
eqns=[diff(Q,x,1)==-(lp(dp-qo*dpio/Q*exp(jjw/(k3*Q^.4))))*am, diff(dp,x,1)==-k2*Q];
conditions=[dp(0)==dpo Q(0)==qo];
dpq=vpa(dsolve(eqns,conditions));
end 


Error using symfun/subsindex (line 157)
Indexing values must be positive integers, logicals,
or symbolic variables.

Error in ro_cp_model (line 4)
eqns=[diff(Q,x,1)==-(lp(dp-qo*dpio/Q*exp(jjw/(k3*Q^.4))))*am,
diff(dp,x,1)==-k2*Q];

Error in mhlro (line 20)
eqns=ro_cp_model(lp,qo,dpio,k3,k2,dpo);
 
>> 

I would like to show a matrix of the anwsers of the two diff. equations with the only uknown value beng jjw.

Answer 1

[Walter Roberson]

lp(dp-qo*dpio/Q*exp(jjw/(k3*Q^.4)))

You have requested variable lp indexed at expression dp-qo*dpio/Q*exp(jjw/(k3*Q^.4)) but that expression involves symbolic variables, and you can never index at a symbolic location.

MATLAB does not have implied multiplication in any context, not even inside the symbolic toolbox.

[Rodrigo Blas]

function dpq=ro_cp_model(lp,qo,dpio,k3,k2,dpo,a1)

syms jjw dp(x) Q(x)
eqns=[diff(Q,x,1)==-(lp*(dp-qo*dpio/Q*exp(jjw/(k3*Q^.4))))*a1, diff(dp,x,1)==-k2*Q];
conditions=[dp(0)==dpo Q(0)==qo];
dpq=vpa(dsolve(eqns,conditions));
end 
clc, clear
a1=28*7*400;                        %%ft^2
po=40.6;                            %%psi
dpio=9;                           %%psi
qf=1525;                            %%gpm
qo=qf;
pp=16.4;                            %%psi
qp=1234;                            %%gpm
depr1=24.9;                         %%psi
delp=depr1-pp;                      %%psi
delpr2=18.2;                        %%psi
PH=360;                                     %%ft
lp=2.76*.264/(60)^2*14.5;                   %%gal/m*min^2/psi
k2=1.75*10^-4*60^1.67*14.5/264.172^1.67;     %%psi/(gal/min)^1.67
k3=2.27*10^-2;                                  %%(m^3/h)^.6/m^2
g=32.17;                            %%ft/s^2
rhoh20=62.4;                        %%lbs/ft^3
pf=rhoh20*g*PH*1.488*.000145;       %%psi
dpo=pf-po;                     %%psi
eqns=ro_cp_model(lp,qo,dpio,k3,k2,dpo,a1);

Warning: Unable to find explicit solution. 
> In dsolve (line 201)
  In ro_cp_model (line 6)
  In mhlro (line 20) 
>> eqns
 
eqns =
 
[ empty sym ]
 
>> 

forgot the multiplication sign. fixed it. Now I get this error

[Walter Roberson]

60^1.67

Really? Or should it be 60^(5/3) ?

.264

So what number is that, exactly ?

• 0.2640000000000000124344978758017532527446746826171875 exactly, which is the closest representable IEEE 754 double precision number ot 0.264 ??
• 1188950301625811/4503599627370496 exactly which is the rational fraction that yields 0.2640000000000000124344978758017532527446746826171875 exactly ?
• 264/1000 exactly, except reduced to lowest terms 33/125 exactly ?
• 264/1000 exactly, but not reduced to lowest terms, because in contexts such as 60^1.67 whether you reduce to lowest terms or not can affect your calculation; for example (-1)^(264/1000) is not the same as (-1)^(33/125)

These things matter, because you are using dsolve(), and dsolve() is for use with calculating exact symbolic solutions if it can. dsolve() is probably going to give up quickly on Q^.4 as that is defined as exp(0.4*log(Q)) which is too messy to dsolve in most cases. But Q^(2/5) is a different matter to the solver, something potentially workable through techniques such as integration by parts, and reducing quintics to quartics by factoring.

If you want to use dsolve() for symbolic solutions, then you should nail down the symbolic representation you want for every number in your equations, converting each of them to rationals or pi or algebraic roots such as sqrt(sym(3)) . If for some reason in your situation it does not make sense to say that 1.488 is exactly 186/125, then you should probably be using a numeric solver instead of dsolve()

[Rodrigo Blas]

Ok I see what you mean. I will try this in ode45. Thanks

[Walter Roberson]

I suggest you look at odeFunction() and in particular the first example which shows a good flow to work through to convert symbolic ODE into something to use numerically.