I'm trying to estimate pi by solving the IVP
$y'' + y = 0$
where
$y(0) = 1, y'(0) = 0$
numerically by defining $\frac{\pi}{2}$ as the first value on t such that $y(t) = 0$
I'm trying to solve this problem both with Euler forward method and with Runge-Kutta 4 method.
However as of right now I get the exact same value from both Euler and RK4 which I think is weird since the error should be higher in the Euler method. (?)
Right now I'm stuck since I don't know what's causing this. I think there could be either one of two problems (maybe both):
- My RK4 is somehow wrong
- The way I try to estimate pi is wrong
I can't find any problem with my RK4 method but I'm not overly confident with this method.. yet.
I think (2) could be the error but if that's the case I don't know how I'm supposed to estimate pi.
I'm aware of the Newton-Raphson method and the bisection method but both of them require a function.. all I have is y values & t values, not a function.
Right now all I'm doing is stopping the simulation as soon as $y(t)$ drops below zero at which I assume t is close enough to $\pi$, which is the first root.
If anyone has any ideas regarding what I'm doing wrong or how I could fix this I'd be very happy.
I'll paste my code (in python) below:
import numpy as np
'''
D.E:
y' = u
u' = -y
'''
def f(u):
return np.array([u[1],-u[0]])
# Euler:
def Euler(y0, t, dt):
t0 = t[0]
for _ in range(1, len(t)):
y0 = y0+dt*f(y0)
if y0[0]<0: # I assume that we're close to pi if y drops below zero and if that's the case we want to exit the loop
return y0,t0
t0 = t0 + dt
return y0, t0
# RK 4
def RK4step(u,dt):
k1 = f(u)
k2 = f(u+0.5*k1*dt)
k3 = f(u+0.5*k2*dt)
k4 = f(u+k3*dt)
return u + dt*(k1+2*k2+2*k3+k4)/6
def RK4(y0, t, dt):
t0 = t[0]
for _ in range(1, len(t)):
y0 = RK4step(y0,dt)
if y0[0] < 0: # I assume that we're close to pi if y drops below zero and if that's the case we want to exit the loop
return y0,t0
t0 = t0 + dt
return y0, t0
if __name__ == "__main__":
real_pi = np.pi
# Initial values
T1 = np.pi/2
T2 = 4
u0 = np.array([1, 0]) # Initial values
n = 5
# Want to loop for different lengths of dt, ideally more than 5 but the program gets very slow after 5.
for i in range(n):
dt = 10**-(i+1)
n = int((T2-T1)/dt)
t = np.linspace(T1,T2,n)
y_euler,t_euler = Euler(u0, t, dt)
y_RK4,t_RK4 = RK4(u0, t, dt)
print()
print(f"{i+1}:")
print(f"------Euler:------")
print(f"Estimated pi: {t_euler}")
print(f"Error Euler = {np.abs(t_euler-real_pi)}")
print(f"------Runge-Kutta 4:------")
print(f"Estimated pi: {t_RK4}")
print(f"Error RK4 = {np.abs(t_RK4-real_pi)}")
```