I have a second order differential equation that I want to solve it in python. The problem is that for one of the variables I don't have the initial condition in 0
but only the value at infinity. Can one tell me what parameters I should provide for scipy.integrate.odeint
? Can it be solved?
Equation:
Theta needs to be found in terms of time. Its first derivative is equal to zero at t=0
. theta is not known at t=0
but it goes to zero at sufficiently large time. all the rest is known. As an approximate I
can be set to zero, thus removing the second order derivative which should make the problem easier.
This is far from being a full answer, but is posted here on the OP's request.
The method I described in the comment is what is known as a shooting method, that allows converting a boundary value problem into an initial value problem. For convenience, I am going to rename your function theta
as y
. To solve your equation numerically, you would first turn it into a first order system, using two auxiliary function, z1 = y
and z2 = y'
, and so your current equation
I y'' + g y' + k y = f(y, t)
would be rewitten as the system
z1' = z2
z2' = f(z1, t) - g z2 - k z1
and your boundary conditions are
z1(inf) = 0
z2(0) = 0
So first we set up the function to compute the derivative of your new vectorial function:
def deriv(z, t) :
return np.array([z[1],
f(z[0], t) - g * z[1] - k * z[0]])
If we had a condition z1[0] = a
we could solve this numerically between t = 0
and t = 1000
, and get the value of y
at the last time as something like
def y_at_inf(a) :
return scipy.integrate.odeint(deriv, np.array([a, 0]),
np.linspace(0, 1000, 10000))[0][-1, 0]
So now all we need to know is what value of a
makes y = 0
at t = 1000
, our poor man's infinity, with
a = scipy.optimize.root(y_at_inf, [1])
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With