GetMeAnIDPls/285-codes
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

fix: RK

Browse Source
This commit is contained in:
fracture-hikari 2025-06-01 21:05:24 +08:00
parent 0b641bdaf7
commit 26d832d64e
1 changed files with 8 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
numeric.py
View File

@ -3,7 +3,7 @@ import matplotlib.pyplot as plt
threshold = 100
def RK_n(y0, t: np.ndarray, n=4):
def RK4(y0, t):
ODE = lambda t, y: t**3 + y**3
#note: t is a vector of time points, y0 is the initial value
steps = len(t)
@ -16,16 +16,12 @@ def RK_n(y0, t: np.ndarray, n=4):
t0_idx = np.where(t >= 0)[0][0]
y[t0_idx] = y0 # set the initial condition at t[0]
#forward RK
for i in range(t0_idx, steps - 1):
#n-th order Runge-Kutta method
k = np.zeros(n)
k[0] = h * ODE(t[i], y[i])
for j in range(1, n):
t_j = t[i] + h * (j / n)
y_j = y[i] + k[j - 1] / n
k[j] = h * ODE(t_j, y_j)
y[i + 1] = y[i] + np.sum(k) / n
k1 = h * ODE(t[i], y[i])
k2 = h * ODE(t[i] + h / 2, y[i] + k1 / 2)
k3 = h * ODE(t[i] + h / 2, y[i] + k2 / 2)
k4 = h * ODE(t[i] + h, y[i] + k3)
y[i + 1] = y[i] + (k1 + 2 * k2 + 2 * k3 + k4) / 6
return y
def euler(f, a, b, y0, h):
@ -84,7 +80,7 @@ for idx, (y0, xm, xM, ym, yM) in enumerate(lims):
Te, Ye = euler(ODE, 0, xM, y0, h)
Ti, Yi = improvedEuler(ODE, 0, xM, y0, h)
Tr = np.linspace(0, xM, int((xM - 0)/h))
Yr = RK_n(y0, Tr)
Yr = RK4(y0, Tr)
axes[idx][0].plot(Te, Ye, label=f'$\Delta t$={h}', linewidth=1)
axes[idx][1].plot(Ti, Yi, label=f'$\Delta t$={h}', linewidth=1)