0
|
6943 Views
|
4 Replies
|
0 Total Likes
View groups...
Share
GROUPS:

# solving systems of equations with given initial conditions

Posted 9 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 == 1.345, nb == 1.412, nw == 0, T == 727.6, P == 2.5}, {na, nb, nw, T, P}, {z, 20}]  is what i thought would work but it is not Attachments:
4 Replies
Sort By:
Posted 9 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 [] that I added, gives In:= Off[NDSolve::ndsz] In:= 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 == 1.345, nb == 1.412, nw == 0, T == 727.6, P == 2.5}, {na, nb, nw, T, P}, {z, 20}][] Out= {na -> InterpolatingFunction[{{0., 2.0480816649058247*^-10}}, <<>>} 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:= naf = na /. sols Out= InterpolatingFunction[{{0., 2.0480816649058247*^-10}},<>] In:= naf[2*10^-10] Out= 1.34499 + 0. I In:= Plot[naf[z], {z, 0, 2*10^-10}] Out= <<>> 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 meNDSolve::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 9 years ago
 here is the modified code... it seems to be something else perhaps? Attachments:
Posted 9 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 includedand eliminate your H[T[z]] :=...Then see what you get.
Posted 9 years ago
 ..