Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
1998
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 1998

[Date Index] [Thread Index] [Author Index]

Search the Archive

RE: DSolve and Erf trouble




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




  • Prev by Date: Re: Re: Formatting posts
  • Next by Date: Random Order of a Set of Numbers
  • Prev by thread: DSolve and Erf trouble
  • Next by thread: Mathematica Palettes