RE: DSolve problems with system of ODEs. Solutions and Thanks
- To: mathgroup at smc.vnet.net
- Subject: [mg21786] RE: DSolve problems with system of ODEs. Solutions and Thanks
- From: "SANCHEZ DE LEON, Guillermo" <gsl at fab.enusa.es>
- Date: Thu, 27 Jan 2000 22:57:09 -0500 (EST)
- Sender: owner-wri-mathgroup at wolfram.com
[My usual E-mail is guillerm at gugu.usal.es]
I posted a few days ago the following question:
> > If I evaluate, in order to obtain the analitycal solution, the System of
> ODEs
> >
> > DSolve[{x'[t] == -0.036 x[t] + 0.0124 y[t] + 0.000035 z[t] + 49.3,
> > y'[t] == 0.0111 x[t] - 0.0286 y[t], z'[t] == 0.0039 x[t] - 0.000035
> > z[t],
> > x[0] == 0, y[0] == 0, z[0] == 0}, {x[t], y[t], z[t]}, t]
> >
> >In version 4.0.0 you receive a massage like this
> > (RootSum::"pfn" :
> > "(DSolve`DSolveDump`dysfunction$6[273147 + <<1>>...
> > is not a pure function.)
In version 2.2 you will receive the correct solutions, and in version 3 you
will receive a error message
I have received several ideas for menber of Mathgroup (Thanks every body, it
a pleasure to have the opportunity of share this group with excelents
Mathematica gurus and belong).
Solution (in Math v4.0.0) forcing the numerical solutions accoding with
Jens-Peer Kuska suggestion
> eqns =
>
> {x'[t] == -0.036 x[t] + 0.0124 y[t] + 0.000035 z[t] + 49.3,
> y'[t] == 0.0111 x[t] - 0.0286 y[t],
> z'[t] == 0.0039 x[t] - 0.000035 z[t], x[0] == 0, y[0] == 0, z[0] == 0};
>
Sol = Chop[Simplify[N[DSolve[eqns], {x[t], y[t], z[t]}, t]]], 10^(-7)]
out[2]:=
{{x[t] -> 1800.097 - 718.4499/E^(0.04466875*t) - 854.97/E^(0.0200*t) -
224.897/E^(0.000030632*t),
y[t] -> 698.6390 + 496.2920/E^(0.04466875163472879*t) -
1108.0984/E^(0.0200356*t) - 87.379/E^(0.000030632*t),
z[t] -> 200582.240249792 + 62.7765894/E^(0.04466875*t) +
166.7142446/E^(0.0200356*t) - 200811.9232/E^(0.0000306322*t)}}
Other solution
> I have deleted the graphic display and most of the output.
>
>
> eqns = {x'[t] == -0.036 x[t] + 0.0124 y[t] + 0.000035 z[t] + 49.3,
> y'[t] == 0.0111 x[t] - 0.0286 y[t],
> z'[t] == 0.0039 x[t] - 0.000035 z[t], x[0] == 0, y[0] == 0, z[0] ==
> 0};
>
> soln = DSolve[ eqns, {x[t], y[t], z[t]}, t];
>
> I am not sure what caused the error message.
> We can get what they indicate by Help Browser > GoTo [RootSum::pfn] with
> Other Information chosen from the category buttons. We get
> "Generated when the first argument in RootSum is not a Function
> expression."
>
> However it looks as if it is serious, since the solution is incorrect:
>
> soln /. t -> 0
>
> {{x[0] -> 1.77331, y[0] -> -0.546333, z[0] -> -0.192132}}
>
> But, in reply to your question about the form of the answer:
>
> RootSum[f, form]represents the sum of form[x]for allx that satisfy the
> polynomial equationf[x] == 0.
> This is very useful notation - for mathematics in general, not just for
> Mathematica.
> In this case we can convert to standard notation (but it takes a while)
>
> ToRadicals[soln]
>
> And, with another wait, to traditional notation
>
> TraditionalForm[%]
>
> Since your input includes inexact numbers, you might accept the muchs
> shorter form.
>
> N[soln]
>
> Or maybe a numerical solution over some range of t (I have used -1 to 1)
>
> Numsoln = NDSolve[eqns, {x[t], y[t], z[t]}, {t, -1, 1}]
>
> This is OK at t=0:
>
> Numsoln /. t -> 0
>
> {{x[0] -> 0., y[0] -> 0., z[0] -> 0.}}
>
> To check on the values throughout the t-range it is convenient to convert
> the solution to one for the functions x,y and z:
>
> FunNumsoln = Nsoln /. f_[t] -> f
>
> and, in eqns change Equal to Subtract:
>
> diffs = eqns /. Equal -> Subtract
>
> Now plot diffs with the functions substitued for x, y, z (with a perfect
> solution each would give constant zero function)
>
> Plot[Evaluate[diffs /. FunNumsoln], {t, -1, 1}, PlotRange -> All];
>
> The maximum absolute error looks to be less than 0.000035.
>
>
> Allan
> ---------------------
> Allan Hayes
> Mathematica Training and Consulting
> Leicester UK
> www.haystack.demon.co.uk
> hay at haystack.demon.co.uk
> Voice: +44 (0)116 271 4198
> Fax: +44 (0)870 164 0565
>
> "