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)}")
```
Set
T1=0, you want to approximate $\pi$ without knowing it. Just multiply the result with 2.This will still give the same result in both methods. The reason is that the angle in the Euler method has a quadratic error, each step has an angle of $\arctan(h)=h+\frac{h^3}3+...$, so the accumulation of $n=t/h$ steps gives $n\arctan(h)=t+\frac{th^2}3+...$ The RK4 step is of course even more exact. Thus both methods pass the zero line after the same number of steps, except in some exceptional cases that were not encountered.
To get a better comparison add an estimate of where in the step the root is. Use a Newton step to refine the intersection point, using that the state vector is $[y,y']$. Using the previous point one could also use the secant root or a reverse cubic interpolation.
Perform the state and time updates at the same time to remain consistent. The loop would then look like
With that RK4 reaches full accuracy of an relative error of $10^{-12}$ in the third step