Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2000
*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 2000

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

Search the Archive

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
-------------------------------------------------------------


  • Prev by Date: cellular automata and OOP
  • Next by Date: Re: How can I teach Mathematica new functions
  • Previous by thread: Re: cellular automata and OOP
  • Next by thread: Re: How can I teach Mathematica new functions