0

I'm trying to solve three ode, one of them is first order and others are second order ode. here are the equations:

my equations

and here is the code:

#note that:
p[0]=y
p[1]=dy/dt
p[2]=q
p[3]=u
p[4]=du/dt

#initial values and requirements :

import numpy as np
import time
from scipy.integrate import odeint

g=1e-25
k=4.14*(10**-5)
t=(10**31)*np.logspace(0.247237,2.35443,10**7)
t=t.tolist()
z0=[2.435*1e26,0,1/3400]

#equations:

def my_eqs1(p,t):
     return[p[1],-2*p[2]*1.4441*(10**-33)*np.sqrt(0.31*(p[2]**-4)*(p[2]+1/3400)+0.69)*p[1]-
       (r*p[2])**(2)*p[0],p[2]*p[2]*1.4441*(10**-33)*np.sqrt(0.31*(p[2]**-4)*(p[2]+1/3400)+0.69)]
       
       

#solving:
start_time=time.perf_counter()
r=1e-29
z1=odeint(my_eqs1,z0,t)
end_time=time.perf_counter()
print(end_time-start_time,"seconds")

In the next step, I have to use p[1] to solve another second order ode, so I tried to redefine my_eqs1 as fallows to solve all of them simultaneously:

def my_eqs2(p,t):
    return[p[1],-2*p[2]*1.4441*(10**-2)*np.sqrt(0.31*(p[2]**-4)*(p[2]+1/3400)+0.69)*p[1]-
           ((1e31*r*p[2])**2)*p[0],p[2]*p[2]*1.4441*(10**-2)*np.sqrt(0.31*(p[2]**-4)*(p[2]+1/3400)+0.69)
           ,1e-62*p[4],-((k**2)+(g*p[1]/k))*p[3]]

and new z0:

z0=[2.435*1e26,0,1/3400,1/np.sqrt(2*k),0.35927359523]

But this time odeint gives this answer which is not good at all.

What is wrong with this code?

By redefining functions and scaling t, (as t->1e31*t) I've got some answers:

t=np.logspace(0.247237,2.35443,10**7)

def my_eqs3(p,t):
        return[p[1],-2*p[2]*1.4441*(10**-2)*np.sqrt(0.31*(p[2]**-4)*(p[2]+1/3400)+0.69)*p[1]-
               ((1e31*r*p[2])**2)*p[0],p[2]*p[2]*1.4441*(10**-2)*np.sqrt(0.31*(p[2]**-4)*(p[2]+1/3400)+0.69)
               ,1e-62*p[4],-((k**2)+(g*p[1]/k))*p[3]]
danial
  • 21
  • 3
  • @LutzLehmann all what I need is to solve the third equation, to do this I have to solve the first and the second one. in codes I tried to solve all of them simultaneously (`my_eqs2`) but got nothing. `my_eqs1` shows what the first and the second equation answers must be. – danial Apr 06 '23 at 07:18
  • Your first system is completely wrong. If the sum of the orders is 5, then the state space has dimension 5, as you enumerated at the top of the code. You can not reduce this to dimension 3 (you could leave out the last equation as the coupling is only downwards there). After that you can run into the problem of scales, in that the default tolerances (mainly the absolute one) become completely inappropriate. – Lutz Lehmann Apr 06 '23 at 07:44
  • @LutzLehmann if by `your first system is completely wrong` you mean the definition of `my_eqs1`, i agree, `my_eqs1` only answers the first and the second equations. if it's wrong, idk how to write it correctly or even what's wrong with it. in `my_eqs2` i just wanted to show how i tried to solve the third equation assuming `my_eqs1` is a correct approach. By changing absolute tolerance i can get answer for the third one, is that what you saying? – danial Apr 06 '23 at 10:55
  • I'm not sure how I mixed them up, I was sure I had seen the third equation in `my_eqc1`. // The best way would be to apply renormalization in time and state to change the components by some factors so that the default tolerances become sensible for the modified system. Otherwise you need to make some trials and build an array of scales and use that in constructing `atol=1e-6*scales`, where as far as I can see `scales = np.array([1e+25,1e-07,1e+00,1e+05,1e+00])` seems to be a sensible initial try. – Lutz Lehmann Apr 06 '23 at 11:12
  • In the second component `dy` you have `,-2*p[2]*1.4441*`, but in the formulas there is no factor `q` at this place. – Lutz Lehmann Apr 06 '23 at 13:12
  • @LutzLehmann thank you, that was a typo, i fixed it. – danial Apr 06 '23 at 21:43

0 Answers0