Home › Forums › Get Help › How do I ...? with Maple

Numerical solution of ode system

psm22's picture
Posted on 2008-07-07 18:27 By psm22 (2)
Categories:

Hi,

I have been trying to get a numerical solution to the following ode system for weeks. I am new to Maple, and I think my issues might be with the coding rather than the maths.

I get the error message:

Error, (in dsolve/numeric/checksing) ode system has a removable singularity at t=0. Initial data is restricted to {v(t) = 0., Phi(t) = 0.}

I know there is a singularity at t=0, so I tried setting the 'range' option to something just above 0, but it is still not working.

Here is my worksheet (without the output, and I am using Maple 10 classic). Sorry, it is pretty ugly:

restart;with(DEtools):with(plots):with(linalg):infolevel[dsolve]:=4:Digits:=20;
> k:=0.065;
> h:=0.5;
> H0:=100*h;
> K:=(omm,omv)->-H0^2*(1-omm-omv);
> Lambda:=(omv)->3*H0^2*omv;
> G:=6.673e-11;#m^3 kg^-1 sec^-2
> a0:=(omm)->25000*omm*h^2:
> a0(1);
> adotarad:=(t)->1/t;
> adotamat:=(t)->(2/t);
t = conformal time:
> arad:=(t)->t/a0(.3):
> amat:=(t)->t^(2)/a0(.3):
> raddom:=loglogplot(arad(t),t=1e-6..1,axes=framed):
> matdom:=loglogplot(amat(t),t=1..a0(.3)^(3/2),axes=framed):
> #display(raddom,matdom);
> rho_dm:=0.3: rho_r:=0.7:

> de1:=diff(Theta0(t),t)=-k*Theta1(t)+(3*adotamat(t))^(-1)*(k^2+3*adotamat(t)^(2)*Phi(t)-4*evalf(Pi)*G*(rho_dm*delta(t)+ 4*rho_r*Theta0(t))):
> de2:=diff(Theta1(t),t)=-k/3*Phi(t)+k/3*Theta0(t):
> de3:=diff(delta(t),t)=-adotamat(t)^(-1)*(k^2+3*adotamat(t))^2*Phi(t)-4*evalf(Pi)*G*rho_dm*delta(t)+ 16*evalf(Pi)*G*rho_dm**rho_r*Theta0(t)-I*k*v(t):
> de4:=diff(v(t),t)=I*k*Phi(t)-adotamat(t)*v(t):
>de5:=diff(Phi(t),t)=3*adotamat(t)*(4*evalf(Pi)*G*amat(t)^2*rho_dm*delta(t)+16*evalf(Pi)*G*amat(t)^2*rho_r*Theta0(t))-(k^2+3*adotamat(t)^2)*Phi(t):

> deqs := de1,de2,de3,de4,de5;
> ics := Theta0(0)=1e-3;Theta1(0)=1e-3;Phi(0)=0.002;v(0)=0.1;delta(0)=0.003;
> #ics := Phi(0)=0.1;Theta0(0)=0.001,Theta1(0)=0.001,delta(0)=0.001,v(0)=0.1;
> fns := Phi(t),Theta0(t),Theta1(t),delta(t),v(t);
> eqsol := dsolve({deqs,ics},{fns},numeric,'range'=1..1e6);

Any help/comments hugely appreciated.

Peter

‹ Runge-Kutta of order 2 does not evaluate 3dplot ›
Doug Meade's picture

Singularity at t=0

Comment on 2008-07-09 10:10 from Doug Meade (Maple Rating 4 721)

Your system of ODEs has a singularity at t=0, so specifying ICs there can be problematic. Are you sure this is the problem you want to solve? You might be able to get away with changing the initial time to a "small" time, but I'll use t=1 just to check that Maple is capable of giving a solution to your system.

While I have not done a full analysis, this system has all the markings of a stiff system. Some coefficients are 1e-10, some are 1e-17, and others are order 1. This is likely to make this problem very difficult to solve.

As for the Maple, the only problem I see is the definition of the initial conditions (ics). You have only one IC, followed by separate equations. Change the semi-colon's to commas, and you will have a full set of ICs.

Making the syntax change and changing the initial time to t=1, Maple will give a numerical solution.

Watching the information from infolevel, I also see that the equations are complex valued? Is this correct?

Here is some output from the numerical solver:

 

eqsol(0.01);
[          
[t = 0.01, 

                                   513                           495    
  Phi(t) = 1.7561641855734908122 10    - 8.8862566109881525003 10    I, 

                                       510                           493    
  Theta0(t) = -2.9293833977001900868 10    + 1.4822789803922300535 10    I, 

                                      506                           489    
  Theta1(t) = 3.1814523334036783304 10    - 1.6098268067697731582 10    I, 

                                     511                           494    
  delta(t) = 2.6364822192224239381 10    - 1.3340698861061906950 10    I, 

                                  489                           506  ]
  v(t) = -4.8295006444616233419 10    - 9.5443969685953777495 10    I]

 

eqsol(0.1);
[                                          44                           26    
[t = 0.1, Phi(t) = 1.6086847331654627340 10   - 8.1400050533538737390 10   I, 

                                       42                           25    
  Theta0(t) = -2.7038561315335253698 10   + 1.3681613382950502379 10   I, 

                                      39                           22    
  Theta1(t) = 3.0043068299750770752 10   - 1.5201905179830742978 10   I, 

                                     43                           26    
  delta(t) = 2.4338163766357651140 10   - 1.2315202100393673506 10   I, 

                                  22                           39  ]
  v(t) = -4.5625739400441621539 10   - 9.0168777453811162121 10   I]

 

eqsol(0.5);
 [                                                                  -15    
 [t = 0.5, Phi(t) = 326.19430771200165313 - 1.6505489550294421852 10    I, 

                                                                -14    
   Theta0(t) = -28.419868131110325123 + 3.4209043642188798897 10    I, 

                                                                -16    
   Theta1(t) = 0.17714921898605595204 - 1.2390596890630632867 10    I, 

   delta(t) = 255.98106937583365131 + 0.0065000000007336814662 I, 

                                                           ]
   v(t) = 0.40000000003012568494 - 0.53612638482822179542 I]

Here's a plot of the real part of each of the 5 components of the solution. Note that I have scaled Theta0 and v to show more detail of the other components. (Call me lazy. I did not re-upload the plot after I changed the labels in the legend.)

odeplot( eqsol, [[t,Re(Phi(t))],[t,Re(Theta0(t))/5],[t,Re(Theta1(t))],[t,Re(delta(t))],[t,Re(v(t))/10]], t=0.5..10, legend=["Phi","Theta0/5","Theta1","delta","v/10"] );

I hope this is helpful.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.ed
psm22's picture

Thanks very much for your

Comment on 2008-07-09 20:56 from psm22 (2)

Thanks very much for your response Doug. The point about trying to start things off in a singularity is clear to me now!

Yes, the solutions should be complex-valued. They represent the evolution of matter density perturbations in the universe, and these oscillated in the early universe when things were much hotter.

I need to think about the time scale issue now.  It is a great help just to have the numerical integration working.

One more quick question: does odeplot have a logarithmic axes option?

Peter

Doug Meade's picture

log plot from odeplot

Comment on 2008-07-09 23:21 from Doug Meade (Maple Rating 4 721)

Peter,

I do not believe odeplot has an explicit log-plot option. But, you can create this by putting the log of whichever component you want to have plotted. One of the nice things about odeplot is that you can specify any combination of dependent and independent variables in the ODE system. I have used this to study conserved quantities, phase plots, .... With a little practice, I'm confident you will be able to produce a plot that suits your needs. If not, post what you have and see if someone on Mapleprimes can help you out.

Good luck,

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.ed

axis=[mode=log]

Comment on 2008-07-10 00:12 from Joe Riel (Maple Rating 4 1258)

You can pass the option axis=[mode=log] to plots[odeplot] and it will use a loglog plot.  Use axis[2]=[mode=log] to affect just the y-axis.  See ?plot,axis

Comment viewing options

Select your preferred way to display the comments and click "Save settings" to activate your changes.
}