Re: Systems of ODEs
- To: mathgroup at smc.vnet.net
- Subject: [mg42082] Re: Systems of ODEs
- From: "Will Self" <wself at msubillings.edu>
- Date: Wed, 18 Jun 2003 02:11:08 -0400 (EDT)
- References: <bcmpdi$spk$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
"PAUL ANDREW ELLISON" <paellison at wisc.edu> wrote in message news:bcmpdi$spk$1 at smc.vnet.net... > Hello everyone, > > I have been solving systems of ODE's with Mathematica recently, and have developed a few questions. To begin with, all the systems of the form A'[t]==C1 A[t]+C2 B[t], B'[t]==C3 A[t]+C4 B[t], and the coefficients (C1, C2, C3, C4) are completely symbolic. The DSolve method worked very well when the problem was simply a system of 2 equations: > > Simplify[DSolve[{P'[t]==k A[t]+k B[t]+k C[t],A'[t]==k A[t],B'[t]==k B[t],C'[t]==k C[t],P[0]==0,A[0]==A0,B[0]==B0,C[0]==C0},{P[t],A[t],B[t],C[t]},t]] > > However, when the system is complicated by adding equations and difficulty to the symbolic coefficients, the calculation takes an extreme amount of time (Calculation was let run for about 70 hours without crashing on a 1.8mHz P4 with 512mb of RAM until i stopped the calculation). The system that caused this long of a calculation was as follows: > > Simplify[DSolve[{P'[t]==k1 A[t]+k2 B[t]+k3 C[t],A'[t]==(-k1-k12) A[t]+k21 B[t],B'[t]==(-k2-k21-k23) B[t]+k12 A[t]+k32 C[t],C'[t]==(-k3-k32) C[t]+k23 B[t],P[0]==0,A[0]==A0,B[0]==B0,C[0]==C0},{P[t],A[t],B[t],C[t]},t]] > > Is there some trick to making such complicated systems solve without having to devote an entire weeks worth of computer time to it? > The only other method that I know to solve them would be to use the eigenvector/eigenvalue form to solve first order systems such as this one. The only problem there is that when I calculate the Eigensystem of the coefficient matrix, the output is nearly 50 pages long. It seems like it would be too much work to copy and paste 50 pages of eigenvectors and eigenvalues into V*e^(\[Lambda]*t) where V=eigenvector and \[Lambda]=corresponding eigenvalue, and then put all the results into the form of the fundamental matrix. Perhaps there is a built-in command for taking the output of the Eigensystem command to a fundamental matrix for the solution of a first order differential system. > Any input would be much appreciated, Thanks in advance. > > -Paul Ellison The fact that the output for the eigensystem was 50 pages long should tell you that you might be expecting too much. Mathematica has to solve a linear system just like everyone else; find the eigenvalues of the coefficient matrix. That means solving a fourth degree equation, and the whole thing will be a mess. I can't imagine being able to use the output for anything, even if you got it. Notice in your case that it is really a 3 by 3 system, because p[t] does not appear on the right side of the equations. So you could solve DSolve[{a'[t] == (-k1 - k12) a[t] + k21 b[t], b'[t] == (-k2 - k21 - k23) b[t] + k12 a[t] + k32 c[t], c'[t] == (-k3 - k32) c[t] + k23 b[t], a[0] == a0, b[0] == b0, c[0] == c0}, {a[t], b[t], c[t]}, t] get that solution (.36 seconds on my computer, Athlon 2200), then use the solutions you got to solve a one-variable equation for p[t]. Everything will be in terms of roots of the characteristic polynomial, though. Notice I've changed your caps to smalls. You are courting trouble to use capital letters, or capitalized symbols, in Mathematica. In fact, look at this: DSolve[C'[t] == C[t], C[t], t] ----> {{C[t] -> E^t*C[1]}} which is just plain wrong. Maybe you can track down why. Hint, it has to do with using capital letters. If something might take a long time, then you would be well advised NOT to wrap it in Simplify. Sometimes Simplify takes a long time. For all you know, your computer had actually found the solution within the 70 hours, and was in the process of running Simplify on the output. Almost always, first get a result, then simplify it as a separate step. You can always abort the Simplify but will still have the result. Presumably you want a symbolic solution so that you can draw general conclusions about the nature of the solutions. I can't imagine any other reason. That strikes me as unlikely because the symbolic solution will be such a mess. Better to see if you can draw any general conclusions about the nature of the eigenvalues, which then you can translate into general conclusions about the solutions. (e.g., if all the eigenvalues have negative real parts, then the solution dies away exponentially with time)