feat: kutta, all numeric and all picard
This commit is contained in:
parent
1a7fd3adf6
commit
0b641bdaf7
116
kutta.py
Normal file
116
kutta.py
Normal file
@ -0,0 +1,116 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
|
||||
def RK_n(y0, t, n=4):
|
||||
ODE = lambda t, y: t**3 + y**3
|
||||
#note: t is a vector of time points, y0 is the initial value
|
||||
steps = len(t)
|
||||
y = np.zeros(steps)
|
||||
h = (t[-1] - t[0]) / (steps - 1)
|
||||
if t[0] > 0:
|
||||
raise ValueError("t[0] must be less than or equal to 0")
|
||||
|
||||
#find the index of the first non-negative time point
|
||||
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
|
||||
|
||||
#backward RK
|
||||
for i in range(t0_idx - 1, -1, -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] = y[i + 1] - np.sum(k) / n
|
||||
return y
|
||||
|
||||
case_num = 3
|
||||
if case_num == 1:
|
||||
plt.title("RK4 result with IVPs for large t")
|
||||
initial_values = [-1, 0, 1]
|
||||
t = np.linspace(0, 2, 1000)
|
||||
for i, y0 in enumerate(initial_values):
|
||||
res = RK_n(y0, t, 4)
|
||||
plt.plot(t, res, label=f'y(0)={y0}, n=4')
|
||||
plt.ylim(-5,5)
|
||||
#observation: the solutions diverge as t increases and seem to have an asymptote which depends on y(0)
|
||||
# and converge to a same solution as t decreases.
|
||||
|
||||
elif case_num == 2:
|
||||
plt.title("RK4 result with IVPs for small t")
|
||||
initial_values = [-1, 0, 1]
|
||||
t = np.linspace(-0.5, 0.5, 1000)
|
||||
for i, y0 in enumerate(initial_values):
|
||||
res = RK_n(y0, t, 4)
|
||||
plt.plot(t, res, label=f'y(0)={y0}, n=4')
|
||||
plt.ylim(-5,5)
|
||||
|
||||
elif case_num == 3:
|
||||
plt.title("Comparing different order RK methods")
|
||||
order = [1, 2, 4, 5]
|
||||
t = np.linspace(-1, 1.4, 1000)
|
||||
res_5 = RK_n(0, t, 5)
|
||||
for n in order:
|
||||
res = RK_n(0, t, n)
|
||||
if n == 1:
|
||||
label = 'Euler - RK5'
|
||||
elif n == 2:
|
||||
label = 'Improved Euler - RK5'
|
||||
else:
|
||||
label = f'RK{n} - RK5'
|
||||
plt.plot(t, res - res_5, label=label)
|
||||
plt.legend()
|
||||
|
||||
elif case_num == 4:
|
||||
plt.title("RK4 with IVP from 0 to 1")
|
||||
initial_values = np.linspace(0, 1, 5)
|
||||
t = np.linspace(-3, 10, 1000)
|
||||
for i, y0 in enumerate(initial_values):
|
||||
res = RK_n(y0, t, 4)
|
||||
plt.plot(t, res, label=f'y(0)={y0}')
|
||||
plt.ylim(-5,10)
|
||||
plt.legend()
|
||||
|
||||
elif case_num == 5:
|
||||
fig, axs = plt.subplots(2, 2, figsize=(10, 8))
|
||||
plt.title("RK4 with different y-scale in the negatives")
|
||||
scale = [5, 10, 20, 50]
|
||||
t = np.linspace(-10, 1, 1000)
|
||||
res = RK_n(0, t, 4)
|
||||
for i, s in enumerate(scale):
|
||||
ax = axs[i // 2, i % 2]
|
||||
ax.plot(t, res / s, label=f'scale={s}')
|
||||
ax.set_ylim(-1, scale[i])
|
||||
ax.axhline(0, color='black', lw=0.5)
|
||||
ax.axvline(0, color='black', lw=0.5)
|
||||
ax.set_title(f'RK4 for t < {s}')
|
||||
# observation: There is no asymptote in the negative y-scale
|
||||
|
||||
elif case_num == 6:
|
||||
plt.title("RK4 with different step sizes")
|
||||
steps = [5, 10, 100, 1000]
|
||||
for n in steps:
|
||||
t = np.linspace(0, 1, n)
|
||||
res = RK_n(0, t, 4)
|
||||
plt.plot(t, res, label=f'{n} steps')
|
||||
plt.legend()
|
||||
|
||||
|
||||
plt.axhline(0, color='black', lw=0.5)
|
||||
plt.axvline(0, color='black', lw=0.5)
|
||||
plt.show()
|
||||
# plt.legend()
|
||||
98
numeric.py
Normal file
98
numeric.py
Normal file
@ -0,0 +1,98 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
threshold = 100
|
||||
|
||||
def RK_n(y0, t: np.ndarray, n=4):
|
||||
ODE = lambda t, y: t**3 + y**3
|
||||
#note: t is a vector of time points, y0 is the initial value
|
||||
steps = len(t)
|
||||
y = np.zeros(steps)
|
||||
h = (t[-1] - t[0]) / (steps - 1)
|
||||
if t[0] > 0:
|
||||
raise ValueError("t[0] must be less than or equal to 0")
|
||||
|
||||
#find the index of the first non-negative time point
|
||||
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
|
||||
return y
|
||||
|
||||
def euler(f, a, b, y0, h):
|
||||
xs = [a]
|
||||
ys = [y0]
|
||||
x = a
|
||||
y = y0
|
||||
while x < b:
|
||||
y += h * f(x, y)
|
||||
x += h
|
||||
if abs(y) > threshold:
|
||||
break
|
||||
xs.append(x)
|
||||
ys.append(y)
|
||||
return xs, ys
|
||||
|
||||
def improvedEuler(f, a, b, y0, h):
|
||||
xs = [a]
|
||||
ys = [y0]
|
||||
|
||||
x = a
|
||||
y = y0
|
||||
while x <= b:
|
||||
slope1 = f(x, y)
|
||||
y_predictor = y + h * slope1
|
||||
slope2 = f(x + h, y_predictor)
|
||||
y = y + (h/2) * (slope1 + slope2)
|
||||
x = x + h
|
||||
if abs(y) > threshold:
|
||||
break
|
||||
xs.append(x)
|
||||
ys.append(y)
|
||||
return xs, ys
|
||||
|
||||
lims = [
|
||||
(0, 1.4, 1.8, 1, threshold),
|
||||
(1, 0.4, 0.6, 1, threshold),
|
||||
(-1, 0.4, 0.6, -threshold, -1),
|
||||
]
|
||||
|
||||
ODE = lambda t, y: t**3 + y**3
|
||||
fig, axes = plt.subplots(3, 3)
|
||||
|
||||
col_titles = ["Euler's Method", "Improved Euler's", "Runge-Kutta"]
|
||||
row_titles = ["y(0) = 0", "y(0) = 1", "y(0) = -1"]
|
||||
|
||||
# Add column titles
|
||||
for col_idx, col_title in enumerate(col_titles):
|
||||
axes[0][col_idx].set_title(col_title, fontsize=12)
|
||||
|
||||
for row_idx, row_title in enumerate(row_titles):
|
||||
axes[row_idx][0].set_ylabel(row_title, fontsize=12, rotation=90, labelpad=20)
|
||||
|
||||
for idx, (y0, xm, xM, ym, yM) in enumerate(lims):
|
||||
for h in [0.1, 0.01, 0.001, 0.0001]:
|
||||
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)
|
||||
|
||||
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)
|
||||
axes[idx][2].plot(Tr, Yr, label=f'$\Delta t$={h}', linewidth=1)
|
||||
for idy in range(3):
|
||||
axes[idx][idy].set_xlim(xm, xM)
|
||||
axes[idx][idy].set_ylim(ym, yM)
|
||||
axes[idx][idy].legend()
|
||||
|
||||
fig.suptitle("Numeric Solutions for Each ODE Near the Asymptote")
|
||||
plt.show()
|
||||
@ -188,35 +188,36 @@ def picardIteration(t0, y0, num_iters = 20):
|
||||
return phi
|
||||
|
||||
lims = [
|
||||
(0, 2, 0, 100),
|
||||
(0, 1, 1, 100),
|
||||
(0, 1, -1, 100)
|
||||
(0, 2, 0, -1, 100),
|
||||
(0, 1, 1, 0, 100),
|
||||
(0, 1, -1, -100, 0)
|
||||
]
|
||||
num_iter = 6
|
||||
|
||||
t0, tM, y0, yM = lims[1]
|
||||
tm, ym = min(t0, tM), min(y0, yM)
|
||||
fig, axes = plt.subplots(2, 3)
|
||||
for idx, (t0, tM, y0, ym, yM) in enumerate(lims):
|
||||
tm = min(t0, tM)
|
||||
picard_result = picardIteration(t0, y0, num_iter)
|
||||
for idp, poly in enumerate(picard_result):
|
||||
np_poly = np.vectorize(poly)
|
||||
ts = np.linspace(t0, tM, 1000)
|
||||
ys = np_poly(ts)
|
||||
axes[0][idx].plot(ts, ys, label=f'#{idp} iteration')
|
||||
axes[0][idx].legend()
|
||||
axes[0][idx].set_xlim(tm, tM)
|
||||
axes[0][idx].set_ylim(ym, yM)
|
||||
|
||||
#last iter coeffs
|
||||
axes[1][idx].scatter(list(range(poly.degree + 1)), poly.coefficients, s=0.5)
|
||||
axes[1][idx].set_yscale('log')
|
||||
|
||||
picard_result = picardIteration(t0, y0, num_iter)
|
||||
'''
|
||||
for idx, poly in enumerate(picard_result):
|
||||
np_poly = np.vectorize(poly)
|
||||
ts = np.linspace(t0, tM, 1000)
|
||||
ys = np_poly(ts)
|
||||
plt.plot(ts, ys)
|
||||
print(idx, poly.degree, poly)
|
||||
'''
|
||||
poly = picard_result[-1]
|
||||
#np_poly = np.vectorize(poly)
|
||||
#ts = np.linspace(t0, tM, 1000)
|
||||
#ys = np_poly(ts)
|
||||
#plt.plot(ts, ys)
|
||||
#print(poly.degree, poly)
|
||||
print(poly.degree)
|
||||
plt.scatter(list(range(poly.degree + 1)), poly.coefficients, s=0.5)
|
||||
plt.yscale('log')
|
||||
#plt.xlim(t0, tM)
|
||||
#plt.ylim(ym, yM)
|
||||
plt.title(f"coefficients of $\\sum c_n t^n$ at y(0) = {y0} of iterations {num_iter}")
|
||||
plt.ylabel("$c_n$")
|
||||
plt.show()
|
||||
col_titles = ["y(0) = 0", "y(0) = 1", "y(0) = -1"]
|
||||
row_titles = ["plots of solutions", "coefficients for $t^n$"]
|
||||
|
||||
# Add column titles
|
||||
for col_idx, col_title in enumerate(col_titles):
|
||||
axes[0][col_idx].set_title(col_title, fontsize=12)
|
||||
|
||||
for row_idx, row_title in enumerate(row_titles):
|
||||
axes[row_idx][0].set_ylabel(row_title, fontsize=12, rotation=90, labelpad=20)
|
||||
plt.show()
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user