# Shooting method for damped oscillator, with Euler method
# based on bisect.cpp & Euler-oscillator.cpp
import numpy as np
M = 0.1 # mass
k = 10 # spring constant
g = 9.8 # gravity
Lambda = 0.7 # damping coefficient
T = 0.1 # final time
dt= 0.001 # time step of Euler method
x0 = 0. # intial position
x1 = 2 # final position
eps=1e-12 # tolerance for machine error
def f(t, x, v):
return -k/M*x +g - Lambda/M*v # particle acceleration
def xfinal(v0): # return x(T) as calculated by Euler method
x=x0
v=v0
for t in np.arange(0, T+eps, dt):
# print(t, x) # output t, x(t)
dx = dt * v # Euler method iteration
dv = dt * f(t, x, v)
x += dx
v += dv
return x
def F(v0): # solving F(v0)=0 (i.e. x(T)=x1) with bisection
return xfinal(v0)-x1
a=0; b=100 # initial range of v0
N=10 # no. of bisection iterations
for i in range(N):
c=(a+b)/2
if F(a)*F(c)<0:
b=c # take left half of the interval
else:
a=c # take right half of the interval
print (a, b)