[Date Index]
[Thread Index]
[Author Index]
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 oneline 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 avalues?

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
 