Re: DSolve problems with system of ODEs. Out in pure notation?
- To: mathgroup at smc.vnet.net
- Subject: [mg21719] Re: [mg21709] DSolve problems with system of ODEs. Out in pure notation?
- From: Bojan Bistrovic <bojanb at physics.odu.edu>
- Date: Sun, 23 Jan 2000 17:58:02 -0500 (EST)
- Sender: owner-wri-mathgroup at wolfram.com
> If I evaluate,
>
> 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]
>
> I receive a massage like this
> (RootSum::"pfn" :
> "(DSolve`DSolveDump`dysfunction$6[273147 + <<1>>...
> is not a pure function.).
> and them a out in pure notation.
>
> But I want to receive the solution in standard notation.
>
> Also the same sentece work fine in Mathematica 2.2.
>
> Can anyone give a hand?
>
Hmmm...
What I get is acctually
Root::"deg":
DSolve`DSolveDump`dysfunction$45 <<1>> has less than 2 root(s) as
a polynomial in DSolve`DSolveDump`dysfunction$45 <<1>>
a few times and then
Power::"infy" : Infinite expression (1/0) encountered.
Just why this DOSEN'T work, I have no idea. What DOES work is NDSolve;
try ploting the solution of
NDSolve[{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,10000}]
in the range {0,300}, {0,500}, {0,1000}, ... You'll notice that the solution
grows pretty fast (exponentially). If you MUST have the analytical solution,
this is what should work: You have the Matrix Equation
r'[t]= A.r[t] + B
with
In[1]:= r[t]={{x[t]},{y[t]},{z[t]}};
In[2]:= A={{-0.036,0.0124,0.000035},{0.0111,-0.0286,0},{0.0039,0,-0.000035}};
In[3]:= B={{49.3},{0},{0}};
In[4]:= r'[t]={{x'[t]},{y'[t]},{z'[t]}};
In[5]:= SetAttributes[Equal,Listable];
In[6]:= r'[t]== A.t[t]+B
Out[6]= {{x'[t]==49.23-0.036 x[t]+0.0124 y[t]+0.000035 z[t]},{
y'[t]==0.0111 x[t]-0.0286 y[t]},{z'[t]==0.0039 x[t]-0.000035 z[t]}}
If you apply the transformation that diagonalizes the matrix A to the whole
equation, you get the three decoupled equations:
In[7]:= S=Transpose[Eigenvectors[A]];
In[8]:= InvS=Inverse[S];
In[10]:= AA=InvX.A.X // Chop
Out[10]= {{-0.0446036,0,0},{0,-0.0200008,0},{0,0,-0.0000306182}}
In[11]:= Eigenvalues[A]
Out[11]= {-0.0446036,-0.0200008,-0.0000306182}
So, what you have next is: InvS.r'[t] = InvS.A.S.InvS.r[t]+InvS.B which yu rewrite as
R'[t]= AA.R[t]+BB
with
In[12]:= R[t] = {{X[t]},{Y[t]},{Z[t]}};
In[13]:= R'[t]= {{X'[t]},{Y'[t]},{Z'[t]}};
In[14]:= BB=InvS.B;
In[15]:= R'[t]= AA.R[t] + BB
Out[14]= {{X'[t]==39.1241-0.0446036 X[t]},
{Y'[t]==-28.3307-0.02000087 Y[t]},
{Z'[t]==6.17106-0.0000306182 Z[t]}}
The initial values remain the same in new coordinates as well; so the new
system is
In[15]:= sols=Union[
DSolve[{X'[t]==39.1241-0.0446036 X[t], X[0] == 0}, {X[t]}, t][[1]],
DSolve[{Y'[t]==-28.3307-0.02000087 Y[t], Y[0] == 0},{Y[t]}, t][[1]],
DSolve[{Z'[t]==6.17106-0.0000306182 Z[t], Z[0] == 0},{Z[t]}, t][[1]]];
This will complain about: PolynomialGCD::"lrgexp";
In the end, you rotate the solutions back:
In[16]:= S.R[t] /. sols //Simplify
Out[16]= {{1806.69 - 718.897 E^(-0.0446036 t) - 861.344 E^(-0.0200009 t) -
226.447 E^(-0.000030618215 t)},
{ 701.193 + 498.622 E^(-0.0446036 t) - 1111.834 E^(-0.0200009 t) -
87.981 E^(-0.0000306182 t)},
{ 201317.448 + 62.9075 E^(-0.0446036 t) + 168.25 E^(-0.02000087 t\) -
201548.605 E^(-0.0000306182 t)}}
When you plot this, it looks the same as the NDSolve result. Judging from the
analytical solution, I'd say that the problem arrises because you have two
different exponentials scales, one of order 10^-2 and the other of 10^-5.
Here's an instructive example: try solving the system
DSolve[{x'[t] == - aa x[t] +bb y[t] + cc z[t] + dd,
y'[t] == ee x[t] - ff y[t], z'[t] ==gg x[t] - cc z[t],
x[0] == 0, y[0] == 0, z[0] == 0}, {x[t], y[t], z[t]}, t]
which is the same as the original one, but the coefficients aren't numerical.
The result I get is:
Eigensystem::"eivec": Unable to find eigenvector for eigenvalue
Root[(bb cc ee - aa cc ff + <<17>>&,1]
a few times and then
Inverse::"sing":
"Matrix {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}} is singular."
It seems that Mathematica tryes to invert the matrix and fails for some
reason. This is in version 3.0.1 undel Linux.
Hope this helps.
Bye, Bojan
-------------------------------------------------------------
Bojan Bistrovic, bojanb at physics.odu.edu
Old Dominion University, Physics Department, Norfolk, VA
-------------------------------------------------------------