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