CODEWITHSUNDEEP
pythonscipyodenumerical-integration
Alex
Alex

Reputation: 437

Scipy ode integrate to unknown limit of t, based on the size of solution

I am modeling charged particles moving through an electromagnetic field and am using scipy ode. The code here is simplified, obviously, but works as an example. The problem I have is that I want to end the integration after a limit on r, not on t. So, integrate dx/dt up to the point where norm(x) > r.

I don't want to just change the function to integrate over r, however, because the position is a function of t. Can I do a definite integral over an unrelated variable or something? 

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def RHS(state, t, Efield, q, mp):

    ds0 = state[3]
    ds1 = state[4]
    ds2 = state[5]
    ds3 = q/mp * Efield*state[0]
    ds4 = q/mp * Efield*state[1]
    ds5 = q/mp * Efield*state[2]

    r = np.linalg.norm((state[0], state[1], state[2]))

    # if r > 30000 then do stop integration.....?

    # return the two state derivatives
    return [ds0, ds1, ds2, ds3, ds4, ds5]


ts = np.arange(0.0, 10.0, 0.1)
state0 = [1.0, 2.0, 3.0, 0.0, 0.0, 0.0]

Efield=1.0
q=1.0
mp=1.0

stateFinal = odeint(RHS, state0, ts, args=(Efield, q, mp))
print(np.linalg.norm(stateFinal[-1,0:2]))

Upvotes: 2

Views: 143

Answers (1)

user6655984
user6655984

Reputation:

You can control the process by performing integration step by step using scipy.integrate.ode. Example: 

from scipy.integrate import ode
t_initial = 0
dt = 0.1
t_max = 10
r = ode(RHS).set_initial_value(state0, t_initial).set_f_params(Efield, q, mp)
solution = [state0]
while r.successful() and r.t < t_max:
    new_val = r.integrate(r.t + dt)
    solution.append(new_val)
    if np.linalg.norm((new_val[0], new_val[1], new_val[2])) > 30000:
        break
print(solution)

Note that the signature of RHS will have to be changed to def RHS(t, state, Efield, q, mp): for ode, the independent variable comes first, unlike in odeint.

The output is the solution computed in the increments of dt of independent variable, up to the time when the loop ends (either because t_max is reached, or integrator fails, or the condition for break was encountered).

Upvotes: 1

Related Questions