If statement inside Do loop
- To: mathgroup at smc.vnet.net
- Subject: [mg4629] If statement inside Do loop
- From: gt1824a at prism.gatech.edu (Heather Mary Hauser)
- Date: Wed, 21 Aug 1996 03:25:43 -0400
- Organization: Georgia Institute of Technology
- Sender: owner-wri-mathgroup at wolfram.com
I apologize for the length of this post, but the problem is driving me crazy.
Mathematica ignores the If statement inside this Do loop. The program runs
fine, no error messages, but it doesn't do what I want it to. Any help on why
this is would be greatly appreciated.
relevant part of program follows:
Do[
sol1 = BesselI[m, z]; (*solution in first region, constant assumed to be 1*)
dsol1 = D[sol1, x];
eq1 = y'[x] + y[x]*(3*DB2/B2 + 1/x - DG/G - (2*V2*m)/(B2*x^2) +
2*(1 + del)*m*(x*B2)^-1) + w[x]*(D2B2/B2 + DB2/B2 - DB2*DG/B2*G
+ (DB2)^2/B2^2 - 2*V2*m*DB2*(x*B2)^-2 - (2*DV2*m)/(B2*x^2) +
(2*V2*m*DG)/(B2*G*x^2) + (2*V2*m)/(B2*x^3) +
(DB2*m*(1 + del))/(x*B2^2) - m^2/x - al^2*G/B2^2);
eq2 = y[x] - w'[x];
bc1 = N[Evaluate[sol1 /. x ->1]]; (*matching pressures across x = 1,
P = (B/k)*w to get initial value for w[x]*)
w'[x] *)
disp1 = -((G*al)/(I*x*B2^2))*((I*x*B1^2)/(al*(4 - B1^2))*dsol1 +
sol1*((2*I*m*B1)/(al*(4 - B1^2)) - (I*B1*x*DB1)/((al*(4 - B1^2))) -
(I/G)*(B1/B2)*(sol1)*((2*V2*m*B2)/(x*al) - B2*x*DB2/al)));
bc2 = N[Evaluate[disp1 /. x-> 1]];
thesol = NDSolve[{eq1 == 0, eq2 == 0, w[1] == bc1, y[1] == bc2},
{w[x], y[x]}, {x, 1, R}];
sol2 = Evaluate[w[x] /. thesol];
dsol2 = Evaluate[y[x] /. thesol];
dsol2p = -(I*al/x)*((I/G)*((-B2^2*x/al)*dsol2 + (2*V2*m*B2*(x*al)^-1 -
B2*x*DB2/al)*sol2) - (((2*V3*m)/(I*B3*x)) -((x*DB3)/(I*B3*al)))*sol2);
bc21 = N[Evaluate[sol2 /. x-> R]];
bc22 = N[Evaluate[dsol2p /. x -> R]];
d2w/dx2 + 1/x dw/dx - (m2/x2 + al2)w*)
Solve for the constant C using expressions for*)
nu = al*x;
sol31 = BesselI[m, nu]; sol32 = BesselK[m, nu];
dsol31 = D[sol1, x]; dsol32 = D[sol32, x];
s31 = N[Evaluate[sol31 /. x -> R]]; s32 = N[Evaluate[sol32 /. x -> R]];
ds31 = N[Evaluate[dsol31 /. x -> R]]; ds32 = N[Evaluate[dsol32 /. x -> R]];
con = N[Evaluate[(ds32*bc21 - s32*bc22/nu)/(s31*ds32 - ds31*s32) /. x -> R]];
recon = N[Evaluate[Re[con]]]; imcon = N[Evaluate[Im[con]]];
lstrl = First[rlist]; lstim = Firts[imlist];
comre = N[Evaluate[lstrl*recon]];
comim = N[Evaluate[lstim*imcon]];
If[comre <= 0 && comim <= 0, Print["instability near", omi],
Print["no"]];
PrependTo[rlist, recon]; PrependTo[imlist, imcon],
{omi, .7107, .7111, .0001}];
Print[TableForm[rlist]]; Print[TableForm[imlist]];
Heather Hauser
Georgia Tech
heather at eas.gatech.edu
==== [MESSAGE SEPARATOR] ====