import scipy.integrate as integrate
import matplotlib.pyplot as plt
import numpy as np
FIG, AX = plt.subplots()
def setup_plot(min_t, max_t, min_y, max_y):
plt.axis([min_t, max_t, min_y, max_y])
AX.set_autoscale_on(False)
plt.grid(True)
def add_plot(t, y, label):
AX.plot(t, y, label=label)
def show_plots():
plt.legend()
plt.show()
def run_solver(solver, initial_values, t_start, t_end, delta_t):
"""
Runs solver for the given interval with step size delta
"""
num_steps = np.floor((t_end-t_start)/delta_t)
t = np.zeros((num_steps, 1))
y = np.zeros((num_steps, 1))
t[0] = t_start
y[0] = initial_values[0]
i = 1
while solver.successful() and i < num_steps:
solver.integrate(solver.t + delta_t)
t[i], y[i] = solver.t, solver.y[0]
i += 1
return t, y
def get_derivatives(a):
"""
Creates array of derivatives representing the ODE
"""
# ODE
# -(1/2)*f'' + (t**4-a)*f=0
# f'' = 2*(t**4-a)*f
# derivatives
# d/dt(y[0]) = y[1]
# d/dt(y[1]) = 2*(t**4-a)*y[0]
derivatives=lambda t, y: np.array([y[1], 2.*(t**4.-a)*y[0]])
return derivatives
def get_solution(a, initial_values, t_start, t_end, delta_t):
"""
Gets numerical solution for the given parameters and interval
"""
# creates LSODA solver
solver = integrate.ode(get_derivatives(a)).set_integrator('lsoda', rtol=1e-12)
solver.set_initial_value(initial_values, t_start)
t, y = run_solver(solver, initial_values, t_start, t_end, delta_t)
return t, y
def update_interval(a_start, a_end, delta_a, initial_values, t_start, t_end, delta_t):
"""
If "a" is smaller than the exact value, the solution blows up to a positive value,
while if "a" is greater than the exact value, solution blows up to a negative value.
The sequence of choices of "a" that causes the solution to blow up positive and
blow up negative defines an interval that contains the exact value of "a".
If we repeatedly call this method with smaller and smaller delta a, we should get
smaller and smaller intervals and therefore better estimates for "a".
"""
a = a_start
t, y = get_solution(a, initial_values, t_start, t_end, delta_t)
previous_a, previous_y = a, y[-1]
while a <= a_end:
a += delta_a
t, y = get_solution(a, initial_values, t_start, t_end, delta_t)
if previous_y/y[-1] < 0:
return previous_a, a, t, y
previous_a, previous_y = a, y[-1]
return a_start, a_end, t, y
def compute_interval(a_start, a_end, delta_a, initial_values, t_start, t_end, delta_t, num_digits):
for _ in range(num_digits):
a_start, a_end, t, y = update_interval(a_start, a_end, delta_a, initial_values, t_start, t_end, delta_t)
label = "%.15f < a < %.15f" % (a_start, a_end)
add_plot(t, y, label)
print(label)
delta_a = delta_a/10.
if __name__ == '__main__':
# set parameters here
t_start, t_end, delta_t = 0., 5, 0.01
num_digits = 6
# setup plot
setup_plot(t_start, t_end, -1, 1)
# experiment with values of "a" here
a = 10
t, y = get_solution(a, [1, 0], t_start, t_end, delta_t)
#add_plot(t, y, "a = %.30f" % a)
# compute interval for a with given accuracy
compute_interval(0.6, 0.7, 0.01, [1, 0], t_start, t_end, delta_t, num_digits) #e0
#compute_interval(2.3, 2.4, 0.01, [0, 1], t_start, t_end, delta_t, num_digits) #e1
#compute_interval(4.6, 4.7, 0.01, [1, 0], t_start, t_end, delta_t, num_digits) #e2
#compute_interval(7.3, 7.4, 0.01, [0, 1], t_start, t_end, delta_t, num_digits) #e3
#compute_interval(10.2, 10.3, 0.01, [1, 0], t_start, t_end, delta_t, num_digits) #e4
#compute_interval(13.3, 13.4, 0.01, [0, 1], t_start, t_end, delta_t, num_digits) #e5
show_plots()
