Message Boards Message Boards

0
|
7118 Views
|
4 Replies
|
0 Total Likes
View groups...
Share
Share this post:

solving systems of equations with given initial conditions

Posted 10 years ago

I would like to solve {na, nb, nw, T, P} at z=20 given these differential equations and initial conditions where these values are given

d := .686;
L := 6.73;
U := 812;
Te := 800;
eP := .6;
OS := .45;
Dp := .0046;
nao := 1.345;
nbo := 1.412;
mwa := .064;
mwb := .032;
mwi := .014;
G := 4*(nao*mwa + nbo*mwb + ni*mwi)/(Pi*d^2);
dens := (4*G*P[z])/(Pi*d*d*8.314*T[z]*(na[z] + nb[z] + nw[z] + ni));
H298 := 2*(-94470) - (2*(-70950));
aa := 5.697;
ab := 6.713;
aw := 12.13;
ba := 0.016;
bb := -0.0000000879;
bw := .00812;
ya := -0.00001185;
yb := .00000417;
yw := 0;
da := .000000003172;
db := -.000000002544;
dw := 0;
H[T[z]] := (H298 ) + (2*aw - 2*aa - 2*ab)*(T[z] - 
     298) + ((2*bw - 2*ba - bb)/2)*(T[z]^2 - 
     298^2) + ((2*yw - 2*ya - yb)/3)*(T[z]^3 - 
     298^3) + ((2*dw - 2*da - db)/4)*(T[z]^4 - 298^4)


NDSolve[{na'[z] == (-2*((Pi*d*d)/4))*r, nb'[z] == (-(Pi*d*d)/4)*r, 
  nw'[z] == 2*((Pi*d*d)/4)*r, 
  T'[z] == 1*((Pi*d*U*(Te - T[z])) - (((Pi*d*d)/4)*r*H))/(na[z]*Cpa + 
       nb[z]*Cpb + nw[z]*Cpw + ni*Cpi), 
  P'[z] == -1*((1 - eP)/
      eP^3)*(G^2/(dens*OS*Dp*9.81))*(((150*(1 - eP)*.32)/(OS*Dp*G)) + 
      1.75), na[0] == 1.345, nb[0] == 1.412, nw[0] == 0, 
  T[0] == 727.6, P[0] == 2.5}, {na, nb, nw, T, P}, {z, 20}]

is what i thought would work but it is not

Attachments:
POSTED BY: dylan keefe
4 Replies
Posted 10 years ago

1: The help page for NDSolve shows NDSolve[{eqns},{vars},{z,min,max}] and you have only {z,num}, but maybe that is ok.

Ignoring that for the moment, evaluating your notebook, note that leading sols= and trailing [[1]] that I added, gives

In[45]:= Off[NDSolve::ndsz]

In[46]:= sols = NDSolve[{na'[z] == (-2*((Pi*d*d)/4))*r, nb'[z] == (-(Pi*d*d)/4)*r, 
    nw'[z] == 2*((Pi*d*d)/4)*r, 
    T'[z] == 1*((Pi*d*U*(Te - T[z])) - (((Pi*d*d)/4)*r*H))/(na[z]*Cpa + nb[z]*Cpb + nw[z]*Cpw + ni*Cpi), 
    P'[z] == -1*((1 - eP)/eP^3)*(G^2/(dens*OS*Dp*9.81))*(((150*(1 - eP)*.32)/(OS*Dp*G)) + 1.75), 
    na[0] == 1.345, nb[0] == 1.412, nw[0] == 0, T[0] == 727.6, 
    P[0] == 2.5}, {na, nb, nw, T, P}, {z, 20}][[1]]

Out[46]= {na -> InterpolatingFunction[{{0.`, 2.0480816649058247`*^-10}},  <<<snip>>>}

The help page for InterpolatingFunction shows Interpolatingfunction[{domain},table] so it appears that your solution has a domain from zero to 2*10^-10. Let's see if we can get results from that.

In[47]:= naf = na /. sols

Out[47]= InterpolatingFunction[{{0.`, 2.0480816649058247`*^-10}},<>]

In[48]:= naf[2*10^-10]

Out[48]= 1.34499 + 0. I

In[49]:= Plot[naf[z], {z, 0, 2*10^-10}]

Out[49]= <<<ApparentlyPerfectlyReasonablePlotSnipped>>>

So it appears that within the domain you are getting results and can plot them.

Even if I correct the syntax to use {z,0,20} inside NDSolve to tell it you want the solution within that domain the result is unchanged.

If I disable your Off[NDSolve::ndsz] so I can see the error messages begin generated by NDSolve it tells me

NDSolve::ndsz: At z == 2.0480816649058247`*^-10, step size is effectively zero; singularity or stiff system suspected. >>

that your system blew up inside NDSolve at 2*10^-10 and exited.

So, can you look at your system, perhaps look at the denominators, check anything else and perhaps see why it would blow up at that point?

Almost always, disabling warnings and errors doesn't make the problem go away, it just hides what is needed to understand them.

POSTED BY: Bill Simpson
Posted 10 years ago

here is the modified code... it seems to be something else perhaps?

Attachments:
POSTED BY: dylan keefe
Posted 10 years ago

NDSolve demands that all variables be assigned values.

r,H,Cpa,Cpb,Cpw,ni,Cpi, have no values.

Unfortunately NDSolve doesn't provide an error message clearly stating undefined variables are included

and eliminate your H[T[z]] :=...

Then see what you get.

POSTED BY: Bill Simpson
Posted 10 years ago

..

POSTED BY: dylan keefe
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract