RE: DSolve and Erf trouble
- To: mathgroup@smc.vnet.net
- Subject: [mg10839] RE: [mg10780] DSolve and Erf trouble
- From: Ersek_Ted%PAX1A@mr.nawcad.navy.mil
- Date: Tue, 10 Feb 1998 21:01:50 -0500
The lines below should give the plot you want.
In[20]:=
sol=Flatten@DSolve[{q'[t]== -a
q[t] r[t], r'[t]== q'[t]+e, q[0]== Q, r[0]== R},{q[t],r[t]},t];
In[21]:=
qa[t_]=q[t]/.sol;
qq[t_]:=qa[SetPrecision[t,17]]
ra[t_]=(r[t]/.sol)/.q[t]->qa[t];
rr[t_]:=ra[SetPrecision[t,17]]
In[22]:=
eq=(q[0]==q[t]/.sol/.t->0);
q[0]=q[0]/.(First@Solve[eq,q[0]]);
a=1/50;Q=50;R=1/1000;e=1;
In[23]:=
Plot[{qq[t],rr[t]},{t,0,100},PlotRange->All]
Out[23]=
-Graphics- (* deleted *)
(1) You need to solve for ( q[0] ). (2) Use exact values for (a, Q, R,
e). (3) Plot insists on using floating point math, so
you have to trick it with SetPrecision[t,17].
Ted Ersek
PS
I don't know about your first question.
|
|two questions in a one-line problem : |
|1/ First part:
|
|An inconsistency came up in a differential equation : |
|sol=Flatten@DSolve[{q'[t]== -a q[t] r[t], r'[t]== -a q[t] r[t]+e,
q[0]== |Q, r[0]== R},{q[t],r[t]},t]
|Part::partw: Part 2 of r'[q] does not exist. Part::partw: Part 2 of
|r'[q] does not exist. Out[90]=
|DSolve[{q'[t] == -(a q[t] r[t]), r'[t] == e - a q[t] r[t], q[0] == Q, |
r[0] == R}, {q[t], r[t]}, t]
|
|... meaning : "I don't work on stuff like this ..." (:-(( ... |what's
r'[q] got to do with it?
|
|but with a little human intervention, sol=Flatten@DSolve[{q'[t]== -a
|q[t] r[t], r'[t]== q'[t]+e, q[0]== Q, r[0]== R},{q[t],r[t]},t] |
|Solve::ifun: Inverse functions are being used by Solve, so some |
solutions may not be found.
|Out[1]=
|{r[t] -> R + e*t - q[0] + q[t],
| q[t] -> (-2*Sqrt[e])/
| (E^((a*t*(e*t + 2*(R - q[0])))/2)* | (Sqrt[e]*(-2/Q +
(Sqrt[a]*E^((a*(R - q[0])^2)/(2*e))*Sqrt[2*Pi]* |
Erf[(Sqrt[a]*(R - q[0]))/(Sqrt[2]*Sqrt[e])])/Sqrt[e]) - |
Sqrt[a]*E^((a*(R - q[0])^2)/(2*e))*Sqrt[2*Pi]* |
Erf[(Sqrt[a]*(R + e*t - q[0]))/(Sqrt[2]*Sqrt[e])]))} |
|... meaning : "got it right this time, that's what I call input" (:-))
|... but what's q[0] doing in there? I told you it was to be
|"Q".
|
|
|2/ Second part.
|
|Now this solution contains a very large Exponential, multiplied with a
|very small difference of two Erf functions. That's asking for trouble
|in Mathematica 3.0.0 since numerical float of Erf[ large ] goes
|haywire. How can it be reduced to a form appropriate for numerical
|evaluation?
|
|a=.02;Q=50;R=.001;e=1.;
|Plot[{qq[t],rr[t]},{t,0,100},PlotRange->All] |
|works ok, but ...
|
|In[104]:=Clear[a,e,R,Q];a=.05;Q=50;R=0.001;e=1.;
|In[105]:=Table[qq[t],{t,0,20}]
|
|Out[105]=
|{50., 594.0556521297816, 6713.817568919678, 72176.73228552992, |
738091.4452270053, 7.179734807316404*^6, 6.643423078153096*^7, |
5.847371602002258*^8, 4.895699508759986*^9, 3.899007663442238*^10, |
2.953783770341961*^11, 2.128573185928179*^12, 1.459095606974182*^13, |
44.06762248565081, 39.04696606830096, 35.83635349893076, |
34.54389468503202, 33.57633898133737, 32.60143669029274, |
31.61976984484918, 30.63888981226333} |
| ... blows up.
|
|how can I explore larger a-values?
|
|wouter.
|Dr. Wouter L. J. MEEUSSEN
|w.meeussen.vdmcc@vandemoortele.be
|eu000949@pophost.eunet.be