Numerical solutions to differential equations

[suyver]

I've got a horrible system of 8 coupled differential equations:

[tex] \frac{\partial}{\partial t} N_0=-R_{0,2}N_0Y_1+\sum_{j=1}^5W_{j,0}N_j-C_{0,1}^{4,2}N_0N_4 [/tex]

[tex] \frac{\partial}{\partial t} N_1=-R_{1,3}N_1Y_1-W_{1,0}N_1+W_{2,1}N_2+W_{4,1}N_4+C_{0,1}^{4,2}N_0N_4 [/tex]

[tex] \frac{\partial}{\partial t} N_2=-R_{2,4}N_2Y_1+R_{0,2}N_0Y_1-W_{2,1}N_2-W_{2,0}N_2+W_{3,2}N_3+C_{0,1}^{4,2}N_0N_4 [/tex]

[tex] \frac{\partial}{\partial t} N_3=-R_{3,5}N_3Y_1+R_{1,3}N_1Y_1-W_{3,2}N_3-W_{3,0}N_3+W_{4,3}N_4 [/tex]

[tex] \frac{\partial}{\partial t} N_4=+R_{2,4}N_2Y_1-W_{4,3}N_4-W_{4,0}N_4+W_{5,4}N_5-C_{0,1}^{4,2}N_0N_4 [/tex]

[tex] \frac{\partial}{\partial t} N_5=+R_{3,5}N_3Y_1-W_{5,4}N_5-W_{5,0}N_5+[/tex]

[tex] \frac{\partial}{\partial t} Y_0=\sum_{j=0}^3R_{j,j+2}N_jY_1+W^{\rm Yb}_{1,0}Y_1[/tex]

[tex] \frac{\partial}{\partial t} Y_1=-\sum_{j=0}^3R_{j,j+2}N_jY_1-W^{\rm Yb}_{1,0}Y_1[/tex]


(easy to copy-paste from the LaTeX-manuscript that I am writing )

Given my boundry conditions (at t=0 all N_i are zero, as is Y_0; Y_1 is some number), I know that I can forget an analytical solution. So, I need a numerical solution for t>0. The only approach that I can think of, is a finite-timestep approach. Using that:

[tex]\frac{\partial f(x)}{\partial x} \equiv \lim_{\Delta x \rightarrow 0}\frac{f(x+\Delta x)-f(x)}{\Delta x}[/tex]

So: I make sure that my timestep [itex]\Delta x[/itex] is small compared to the (expected) fluctuations in the solutions. In that case, I can re-write the equations in an iterative form (note that I used [itex]\Delta x=1[/itex]):

[tex] f(n+1) = {\rm some\; function\; of}\; f(n)[/tex]

which I can do for all of my functions. Now, starting with the boundry values (for time [itex]n=0[/itex]) I can generate the solutions for [itex]n>0[/itex] with a computer. Computationally intensive, but that should not be a major problem...

However, are there better / faster ways to do this? Since I'd need to do it quite a number of times, a significant time-saver would be appreciated...



Physical background: I am modelling energy-transfer processes in rare-earth doped insulators. The N_i are populations of specific energy levels in the rare Earth ion and the Y_i are populations of a sensitizer ion. Very cool stuff (I have lots of experimental data that I try to put into a "simple" model; see above ).

 

[Dr Transport]

Suyver,

What you have is a mess, nonlinearity abounds. I suggest that you do a stability analysis on this model. Try solving the linear problem first, I would suspect that you will get a solution of the form [tex] e^{-\lambda} [/tex] where [tex] \lambda [/tex] is an eigenvalue of the "linear" matrix. I would go to the library and get a book on mathematical modeling, I looked at my copy of Haberman and found a nonlinear set of equations for population dynamics, which is essentially what this problem is. Add in the nonlinear terms one at a time and see where they lead you.

This seems like a very interestng problem and I would like to see the outcome. Let us know how it is going either on this forum or to personal messages here. I have a couple of advanced degrees in Physics, a Semiconductor modeling PhD, so this intrigues me to the point of trying do a little programming myself, just to see how the results evolve.

Dr Transport

 

[suyver]

Hi, thanks for the response.

I don't get the use of trying to find analytical solutions to a part of the problem first. It seems to be a lot of work and I don't see what use they'd have. After all, I am only interested in the numerical outcome (I want to fit the results, to obtain values for some of the [itex]R_{i,j}[/itex], [itex]W_{i,j}[/itex] as well as for [itex]C_{0,1}^{4,2}[/itex]). I basically already 'know' what the solutions for the N_i look like graphically: they start at 0, then there is a rise to a maximum at a time that I roughly know, and then there is a decay back to zero.

I'll look in the library, but we are rather poorly stocked on numerical mathematics. I'd really hoped that there was an easy numerical method (like the 'finite timestep' method that I described in my first post) that I could use...

Cheers!

PS: PhD in solid state physics here too!

 

[Haelfix]

Suyver, I suggest using Matlab for this problem, they have several different numerical integration packages, one or two might work.

Naively, I would have done exactly what you did, that is the finite-step method.

Runge-Kutta will stall badly if you want any precision, but still computers are pretty powerful these days, and it is of course the easiest to setup.

The reason you want an analytic solution for some subset of the solution space is to check the stability of your system. Its worth it, believe me.

I've been confronted in the past, with systems that output numerical answers that were essentially complete garbage, even with very advanced numerical tools.