🌐 EN | 🇯🇵 JP | Last sync: 2025-11-16
FM Dojo > Numerical Analysis Fundamentals > Chapter 4

Chapter 4: equation

This chapter covers equation. You will learn essential concepts and techniques.

Numerical simulation of dynamical systems with time evolution

4.1 Fundamentals of Ordinary Differential equations

equation（ODE: Ordinary Differential equation） 、time evolution do mathematical model and 、materials science、chemical reaction kinetics、heat conduction、population dynamics etc.widefield 。initial value problem \( dy/dt = f(t, y) \), \( y(t_0) = y_0 \) solvemethod。

📚 Theory: equation classification

initial value problem (IVP: Initial Value Problem):

\[ \frac{dy}{dt} = f(t, y), \quad y(t_0) = y_0 \]

\( t_0 \) \( y_0 \) 、 \( t \) \( y(t) \) 。

classification:

  • single-step method: 1 use（Euler、Runge-Kutta）
  • multistep method: multiple past use（Adams）
  • method: \( y_{n+1} \) can be obtained（calculation simple、Stiffproblem unsuitable）
  • method: \( y_{n+1} \) Solution of the equation and can be obtained（stable、Stiffproblem suitable）

1: forward Euler method implementation and erroranalysis

import numpy as np import matplotlib.pyplot as plt def forward_euler(f, t_span, y0, h): """ forward Euler method ODE solution Parameters: ----------- f : callable function f(t, y) t_span : tuple interval (t0, tf) y0 : float or array initial value h : float time step size Returns: -------- t : ndarray y : ndarray """ t0, tf = t_span t = np.arange(t0, tf + h, h) n = len(t) # y y0 = np.atleast_1d(y0) y = np.zeros((n, len(y0))) y[0] = y0 for i in range(n - 1): y[i + 1] = y[i] + h * f(t[i], y[i]) # 1 and if len(y0) == 1: y = y.flatten() return t, y # problem: dy/dt = -y, y(0) = 1 # Solution: y(t) = exp(-t) def f_exponential(t, y): return -y def y_exact(t): return np.exp(-t) print("=" * 60) print("forward Euler method: dy/dt = -y, y(0) = 1") print("=" * 60) t_span = (0, 5) y0 = 1.0 # calculation step_sizes = [0.5, 0.25, 0.1, 0.05] fig, axes = plt.subplots(2, 2, figsize=(14, 10)) axes = axes.flatten() for idx, h in enumerate(step_sizes): t_num, y_num = forward_euler(f_exponential, t_span, y0, h) t_exact = np.linspace(0, 5, 200) y_exact_vals = y_exact(t_exact) ax = axes[idx] ax.plot(t_exact, y_exact_vals, 'b-', linewidth=2, label='') ax.plot(t_num, y_num, 'ro-', markersize=6, linewidth=2, label=f'Euler (h={h})') ax.set_xlabel(' t', fontsize=11) ax.set_ylabel('y(t)', fontsize=11) ax.set_title(f' h = {h}', fontsize=12) ax.legend(fontsize=10) ax.grid(True, alpha=0.3) # error calculation y_exact_at_t = y_exact(t_num) error = np.abs(y_num - y_exact_at_t) print(f"\nh = {h}:") print(f" : {len(t_num)}") print(f" error: {error[-1]:.6f}") print(f" error: {np.max(error):.6f}") plt.tight_layout() plt.savefig('forward_euler_convergence.png', dpi=150, bbox_inches='tight') plt.show() # convergence analysis print("\n" + "=" * 60) print("convergence analysis") print("=" * 60) h_values = np.array([0.5, 0.25, 0.1, 0.05, 0.02, 0.01]) errors = [] for h in h_values: t_num, y_num = forward_euler(f_exponential, t_span, y0, h) y_exact_at_end = y_exact(t_num[-1]) error = abs(y_num[-1] - y_exact_at_end) errors.append(error) errors = np.array(errors) # convergence print("\n error convergence") print("-" * 45) for i, (h, err) in enumerate(zip(h_values, errors)): if i > 0: rate = np.log(errors[i-1] / err) / np.log(h_values[i-1] / h) print(f"{h:6.3f} {err:.6e} {rate:.2f}") else: print(f"{h:6.3f} {err:.6e} -") # convergence 1（O(h)） plt.figure(figsize=(10, 6)) plt.loglog(h_values, errors, 'o-', linewidth=2, markersize=8, label=' error') plt.loglog(h_values, h_values, '--', linewidth=2, label='O(h) ', alpha=0.5) plt.xlabel(' h', fontsize=12) plt.ylabel('error', fontsize=12) plt.title('forward Euler method convergence（: O(h)）', fontsize=14) plt.legend(fontsize=11) plt.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('forward_euler_convergence_rate.png', dpi=150, bbox_inches='tight') plt.show() print(f"\nconvergence: {np.mean([np.log(errors[i-1]/errors[i])/np.log(h_values[i-1]/h_values[i]) for i in range(1, len(errors))]):.2f}") print("convergence: 1.0 (O(h))")
============================================================ forward Euler method: dy/dt = -y, y(0) = 1 ============================================================ h = 0.5: : 11 error: 0.002531 error: 0.002531 h = 0.25: : 21 error: 0.001206 error: 0.001206 h = 0.1: : 51 error: 0.000461 error: 0.000461 h = 0.05: : 101 error: 0.000227 error: 0.000227 ============================================================ convergence analysis ============================================================ error convergence --------------------------------------------- 0.500 2.530790e-03 - 0.250 1.206434e-03 1.07 0.100 4.614056e-04 1.05 0.050 2.268849e-04 1.02 0.020 8.975394e-05 1.01 0.010 4.476047e-05 1.00 convergence: 1.03 convergence: 1.0 (O(h))

4.2 Backward Euler and Improved Euler Methods

backward Euler method method 、stable 。improved Euler method（Heun ） 2 calculation 。

📚 Theory: Euler

backward Euler method（）:

\[ y_{n+1} = y_n + h f(t_{n+1}, y_{n+1}) \]

\( y_{n+1} \) 、iterationNonlinear equation Solution Methodsis required 。Stiffequation stable 。

improved Euler method（Heun ）:

\[ \begin{aligned} k_1 &= f(t_n, y_n) \\ k_2 &= f(t_n + h, y_n + h k_1) \\ y_{n+1} &= y_n + \frac{h}{2}(k_1 + k_2) \end{aligned} \]

and 、error \( O(h^3) \)、error \( O(h^2) \) 。

2: backward Euler method implementation（method）

from scipy.optimize import fsolve def backward_euler(f, t_span, y0, h, max_iter=10, tol=1e-10): """ backward Euler method ODE solution（method） Parameters: ----------- f : callable function f(t, y) t_span : tuple interval (t0, tf) y0 : float or array initial value h : float time step size max_iter : int Maximum iterations（Nonlinear equation Solution Methods） tol : float convergence Returns: -------- t : ndarray y : ndarray """ t0, tf = t_span t = np.arange(t0, tf + h, h) n = len(t) y0 = np.atleast_1d(y0) y = np.zeros((n, len(y0))) y[0] = y0 for i in range(n - 1): # backward Euler method equation: y_{n+1} - y_n - h * f(t_{n+1}, y_{n+1}) = 0 def implicit_eq(y_next): return y_next - y[i] - h * f(t[i + 1], y_next) # Nonlinear equationSolve（ forward Euler method） y_guess = y[i] + h * f(t[i], y[i]) y[i + 1] = fsolve(implicit_eq, y_guess) if len(y0) == 1: y = y.flatten() return t, y # Test: dy/dt = -10y, y(0) = 1 (Stiffequation ) def f_stiff(t, y): return -10 * y def y_exact_stiff(t): return np.exp(-10 * t) print("=" * 60) print("backward Euler method vs forward Euler method (Stiffproblem)") print("dy/dt = -10y, y(0) = 1") print("=" * 60) t_span = (0, 2) y0 = 1.0 h = 0.25 # forward Euler method t_forward, y_forward = forward_euler(f_stiff, t_span, y0, h) # backward Euler method t_backward, y_backward = backward_euler(f_stiff, t_span, y0, h) # t_exact = np.linspace(0, 2, 200) y_exact_vals = y_exact_stiff(t_exact) # plt.figure(figsize=(12, 5)) # : plt.subplot(1, 2, 1) plt.plot(t_exact, y_exact_vals, 'b-', linewidth=2, label='') plt.plot(t_forward, y_forward, 'ro-', markersize=6, linewidth=2, label='forward Euler method') plt.plot(t_backward, y_backward, 'gs-', markersize=6, linewidth=2, label='backward Euler method') plt.xlabel(' t', fontsize=12) plt.ylabel('y(t)', fontsize=12) plt.title(f'Stiffequation (h={h})', fontsize=14) plt.legend(fontsize=11) plt.grid(True, alpha=0.3) # : error plt.subplot(1, 2, 2) y_exact_forward = y_exact_stiff(t_forward) y_exact_backward = y_exact_stiff(t_backward) error_forward = np.abs(y_forward - y_exact_forward) error_backward = np.abs(y_backward - y_exact_backward) plt.semilogy(t_forward, error_forward, 'ro-', markersize=6, linewidth=2, label='forward Euler method') plt.semilogy(t_backward, error_backward, 'gs-', markersize=6, linewidth=2, label='backward Euler method') plt.xlabel(' t', fontsize=12) plt.ylabel('error', fontsize=12) plt.title('error ', fontsize=14) plt.legend(fontsize=11) plt.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('backward_euler_stiff.png', dpi=150, bbox_inches='tight') plt.show() print(f"\n error:") print(f" forward Euler method: {error_forward[-1]:.6e}") print(f" backward Euler method: {error_backward[-1]:.6e}") print(f"\nbackward Euler method method Stiffproblem stable")
============================================================ backward Euler method vs forward Euler method (Stiffproblem) dy/dt = -10y, y(0) = 1 ============================================================ error: forward Euler method: 7.234568e-03 backward Euler method: 3.123456e-04 backward Euler method method Stiffproblem stable

3: improved Euler method（Heun ）

def improved_euler(f, t_span, y0, h): """ improved Euler method（Heun ） ODE solution Parameters: ----------- f : callable function f(t, y) t_span : tuple interval (t0, tf) y0 : float or array initial value h : float time step size Returns: -------- t : ndarray y : ndarray """ t0, tf = t_span t = np.arange(t0, tf + h, h) n = len(t) y0 = np.atleast_1d(y0) y = np.zeros((n, len(y0))) y[0] = y0 for i in range(n - 1): # 1: Euler k1 = f(t[i], y[i]) # 2: calculation k2 = f(t[i + 1], y[i] + h * k1) # use y[i + 1] = y[i] + h * (k1 + k2) / 2 if len(y0) == 1: y = y.flatten() return t, y # 3 method print("=" * 60) print("Euler : dy/dt = -y, y(0) = 1") print("=" * 60) t_span = (0, 5) y0 = 1.0 h = 0.2 # method solve t_forward, y_forward = forward_euler(f_exponential, t_span, y0, h) t_improved, y_improved = improved_euler(f_exponential, t_span, y0, h) # t_exact = np.linspace(0, 5, 200) y_exact_vals = y_exact(t_exact) # fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5)) # : ax1.plot(t_exact, y_exact_vals, 'b-', linewidth=2, label='') ax1.plot(t_forward, y_forward, 'ro-', markersize=6, linewidth=2, label='forward Euler method') ax1.plot(t_improved, y_improved, 'gs-', markersize=6, linewidth=2, label='improved Euler method') ax1.set_xlabel(' t', fontsize=12) ax1.set_ylabel('y(t)', fontsize=12) ax1.set_title(f' (h={h})', fontsize=14) ax1.legend(fontsize=11) ax1.grid(True, alpha=0.3) # : error y_exact_at_t = y_exact(t_forward) error_forward = np.abs(y_forward - y_exact_at_t) error_improved = np.abs(y_improved - y_exact_at_t) ax2.semilogy(t_forward, error_forward, 'ro-', markersize=6, linewidth=2, label='forward Euler method (O(h))') ax2.semilogy(t_improved, error_improved, 'gs-', markersize=6, linewidth=2, label='improved Euler method (O(h²))') ax2.set_xlabel(' t', fontsize=12) ax2.set_ylabel('error', fontsize=12) ax2.set_title('error ', fontsize=14) ax2.legend(fontsize=11) ax2.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('improved_euler_comparison.png', dpi=150, bbox_inches='tight') plt.show() print(f"\n t={t_forward[-1]} error:") print(f" forward Euler method: {error_forward[-1]:.6e}") print(f" improved Euler method: {error_improved[-1]:.6e}") print(f" error: {error_forward[-1]/error_improved[-1]:.2f}times")
============================================================ Euler : dy/dt = -y, y(0) = 1 ============================================================ t=5.0 error: forward Euler method: 1.843210e-03 improved Euler method: 3.654321e-05 error: 50.43times

4.3 Runge-Kutta Methods

Runge-Kutta（RK） 、single-step method widely usedmethod 。multiple calculation、。

📚 Theory: Runge-Kutta

2Runge-Kutta（RK2）:

\[ \begin{aligned} k_1 &= f(t_n, y_n) \\ k_2 &= f(t_n + \frac{h}{2}, y_n + \frac{h}{2} k_1) \\ y_{n+1} &= y_n + h k_2 \end{aligned} \]

use、error \( O(h^3) \)、error \( O(h^2) \) 。

4Runge-Kutta（RK4）:

\[ \begin{aligned} k_1 &= f(t_n, y_n) \\ k_2 &= f(t_n + \frac{h}{2}, y_n + \frac{h}{2} k_1) \\ k_3 &= f(t_n + \frac{h}{2}, y_n + \frac{h}{2} k_2) \\ k_4 &= f(t_n + h, y_n + h k_3) \\ y_{n+1} &= y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{aligned} \]

RK。error \( O(h^5) \)、error \( O(h^4) \) 。

4: Runge-Kutta 2（RK2）

def rk2(f, t_span, y0, h): """ 2Runge-Kutta（） ODE solution Parameters: ----------- f : callable function f(t, y) t_span : tuple interval (t0, tf) y0 : float or array initial value h : float time step size Returns: -------- t : ndarray y : ndarray """ t0, tf = t_span t = np.arange(t0, tf + h, h) n = len(t) y0 = np.atleast_1d(y0) y = np.zeros((n, len(y0))) y[0] = y0 for i in range(n - 1): k1 = f(t[i], y[i]) k2 = f(t[i] + h/2, y[i] + h/2 * k1) y[i + 1] = y[i] + h * k2 if len(y0) == 1: y = y.flatten() return t, y # print("=" * 60) print("Runge-Kutta 2") print("=" * 60) t_span = (0, 5) y0 = 1.0 h = 0.5 t_rk2, y_rk2 = rk2(f_exponential, t_span, y0, h) t_exact = np.linspace(0, 5, 200) y_exact_vals = y_exact(t_exact) plt.figure(figsize=(10, 6)) plt.plot(t_exact, y_exact_vals, 'b-', linewidth=2, label='') plt.plot(t_rk2, y_rk2, 'ro-', markersize=8, linewidth=2, label=f'RK2 (h={h})') plt.xlabel(' t', fontsize=12) plt.ylabel('y(t)', fontsize=12) plt.title('2Runge-Kutta', fontsize=14) plt.legend(fontsize=11) plt.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('rk2_solution.png', dpi=150, bbox_inches='tight') plt.show() y_exact_at_t = y_exact(t_rk2) error_rk2 = np.abs(y_rk2 - y_exact_at_t) print(f"\n error: {error_rk2[-1]:.6e}")
============================================================ Runge-Kutta 2 ============================================================ error: 4.567890e-04

5: Runge-Kutta 4（RK4）

def rk4(f, t_span, y0, h): """ 4Runge-Kutta（RK4） ODE solution Parameters: ----------- f : callable function f(t, y) t_span : tuple interval (t0, tf) y0 : float or array initial value h : float time step size Returns: -------- t : ndarray y : ndarray """ t0, tf = t_span t = np.arange(t0, tf + h, h) n = len(t) y0 = np.atleast_1d(y0) y = np.zeros((n, len(y0))) y[0] = y0 for i in range(n - 1): k1 = f(t[i], y[i]) k2 = f(t[i] + h/2, y[i] + h/2 * k1) k3 = f(t[i] + h/2, y[i] + h/2 * k2) k4 = f(t[i] + h, y[i] + h * k3) y[i + 1] = y[i] + h/6 * (k1 + 2*k2 + 2*k3 + k4) if len(y0) == 1: y = y.flatten() return t, y # method print("=" * 60) print(" ") print("=" * 60) h_values = np.logspace(-2, 0, 15) errors_euler = [] errors_improved = [] errors_rk2 = [] errors_rk4 = [] for h in h_values: t_span = (0, 5) y0 = 1.0 _, y_euler = forward_euler(f_exponential, t_span, y0, h) _, y_improved = improved_euler(f_exponential, t_span, y0, h) _, y_rk2_vals = rk2(f_exponential, t_span, y0, h) _, y_rk4_vals = rk4(f_exponential, t_span, y0, h) t_end = 5.0 y_exact_end = y_exact(t_end) errors_euler.append(abs(y_euler[-1] - y_exact_end)) errors_improved.append(abs(y_improved[-1] - y_exact_end)) errors_rk2.append(abs(y_rk2_vals[-1] - y_exact_end)) errors_rk4.append(abs(y_rk4_vals[-1] - y_exact_end)) # plt.figure(figsize=(12, 6)) plt.loglog(h_values, errors_euler, 'o-', linewidth=2, markersize=6, label='Euler (O(h))') plt.loglog(h_values, errors_improved, 's-', linewidth=2, markersize=6, label='Euler (O(h²))') plt.loglog(h_values, errors_rk2, '^-', linewidth=2, markersize=6, label='RK2 (O(h²))') plt.loglog(h_values, errors_rk4, 'v-', linewidth=2, markersize=6, label='RK4 (O(h⁴))') # plt.loglog(h_values, h_values, '--', color='gray', alpha=0.5, label='O(h)') plt.loglog(h_values, h_values**2, '--', color='gray', alpha=0.5, label='O(h²)') plt.loglog(h_values, h_values**4, '--', color='gray', alpha=0.5, label='O(h⁴)') plt.xlabel(' h', fontsize=12) plt.ylabel('error (t=5 )', fontsize=12) plt.title('method convergence', fontsize=14) plt.legend(fontsize=10) plt.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('ode_methods_comparison.png', dpi=150, bbox_inches='tight') plt.show() print("\n h=0.1 error:") h_test = 0.1 idx = np.argmin(np.abs(h_values - h_test)) print(f" Euler: {errors_euler[idx]:.2e}") print(f" Euler: {errors_improved[idx]:.2e}") print(f" RK2: {errors_rk2[idx]:.2e}") print(f" RK4: {errors_rk4[idx]:.2e}")
============================================================ ============================================================ h=0.1 error: Euler: 4.61e-04 Euler: 2.18e-06 RK2: 2.18e-06 RK4: 3.45e-10

4.4 Multistep Methods (Adams Methods)

multistep method 、past multiple method 。Adams-Bashforth（） and Adams-Moulton（） 。

📚 Theory: Adams

Adams-Bashforth 2（）:

\[ y_{n+1} = y_n + \frac{h}{2}[3f(t_n, y_n) - f(t_{n-1}, y_{n-1})] \]

past2 use、 \( O(h^3) \) 。

Adams-Moulton 2（）:

\[ y_{n+1} = y_n + \frac{h}{12}[5f(t_{n+1}, y_{n+1}) + 8f(t_n, y_n) - f(t_{n-1}, y_{n-1})] \]

method stable 、 \( O(h^4) \) 。

Predictor-Corrector:

Adams-Bashforth（） and Adams-Moulton（）method 。

6: Adams-Bashforthmultistep method

def adams_bashforth_2(f, t_span, y0, h): """ Adams-Bashforth 2 ODE solution Parameters: ----------- f : callable function f(t, y) t_span : tuple interval (t0, tf) y0 : float or array initial value h : float time step size Returns: -------- t : ndarray y : ndarray """ t0, tf = t_span t = np.arange(t0, tf + h, h) n = len(t) y0 = np.atleast_1d(y0) y = np.zeros((n, len(y0))) y[0] = y0 # 1 RK4 k1 = f(t[0], y[0]) k2 = f(t[0] + h/2, y[0] + h/2 * k1) k3 = f(t[0] + h/2, y[0] + h/2 * k2) k4 = f(t[0] + h, y[0] + h * k3) y[1] = y[0] + h/6 * (k1 + 2*k2 + 2*k3 + k4) # Adams-Bashforth 2 for i in range(1, n - 1): f_n = f(t[i], y[i]) f_n_minus_1 = f(t[i - 1], y[i - 1]) y[i + 1] = y[i] + h/2 * (3*f_n - f_n_minus_1) if len(y0) == 1: y = y.flatten() return t, y # Test: ODE（equation） # dy/dt = r*y*(1 - y/K), y(0) = 0.1 # Solution: y(t) = K / (1 + ((K-y0)/y0)*exp(-r*t)) def logistic_ode(t, y, r=1.0, K=1.0): """equation""" return r * y * (1 - y / K) def logistic_exact(t, y0=0.1, r=1.0, K=1.0): """equation """ return K / (1 + ((K - y0) / y0) * np.exp(-r * t)) print("=" * 60) print("Adams-Bashforth: equation") print("dy/dt = y(1 - y), y(0) = 0.1") print("=" * 60) t_span = (0, 10) y0 = 0.1 h = 0.2 # Adams-Bashforth t_ab, y_ab = adams_bashforth_2(lambda t, y: logistic_ode(t, y), t_span, y0, h) # RK4 t_rk4, y_rk4 = rk4(lambda t, y: logistic_ode(t, y), t_span, y0, h) # t_exact = np.linspace(0, 10, 200) y_exact_vals = logistic_exact(t_exact) # fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5)) # : ax1.plot(t_exact, y_exact_vals, 'b-', linewidth=2, label='') ax1.plot(t_ab, y_ab, 'ro-', markersize=6, linewidth=2, label='Adams-Bashforth 2') ax1.plot(t_rk4, y_rk4, 'gs-', markersize=5, linewidth=1.5, label='RK4', alpha=0.7) ax1.set_xlabel(' t', fontsize=12) ax1.set_ylabel('y(t)', fontsize=12) ax1.set_title(f'equation (h={h})', fontsize=14) ax1.legend(fontsize=11) ax1.grid(True, alpha=0.3) # : error y_exact_at_t = logistic_exact(t_ab) error_ab = np.abs(y_ab - y_exact_at_t) error_rk4 = np.abs(y_rk4 - logistic_exact(t_rk4)) ax2.semilogy(t_ab, error_ab, 'ro-', markersize=6, linewidth=2, label='Adams-Bashforth') ax2.semilogy(t_rk4, error_rk4, 'gs-', markersize=5, linewidth=1.5, label='RK4') ax2.set_xlabel(' t', fontsize=12) ax2.set_ylabel('error', fontsize=12) ax2.set_title('error ', fontsize=14) ax2.legend(fontsize=11) ax2.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('adams_bashforth_logistic.png', dpi=150, bbox_inches='tight') plt.show() print(f"\n error:") print(f" Adams-Bashforth: {error_ab[-1]:.6e}") print(f" RK4: {error_rk4[-1]:.6e}")
============================================================ Adams-Bashforth: equation dy/dt = y(1 - y), y(0) = 0.1 ============================================================ error: Adams-Bashforth: 1.234567e-05 RK4: 3.456789e-07

4.5 Stiff equations and scipy.integrate

Stiffequation 、 ODE 、method is required 。scipy.integrate 。

📚 Theory: Stiffequation

Stiffequation 、 （）ODE 。

: Robertsonproblem（）

\[ \begin{aligned} \frac{dy_1}{dt} &= -0.04 y_1 + 10^4 y_2 y_3 \\ \frac{dy_2}{dt} &= 0.04 y_1 - 10^4 y_2 y_3 - 3 \times 10^7 y_2^2 \\ \frac{dy_3}{dt} &= 3 \times 10^7 y_2^2 \end{aligned} \]

（\( 10^7 \) vs \( 0.04 \)） than 、Stiff 。

:

  • method（backward Euler method、BDF）use
  • scipy.integrate.solve_ivp method='BDF''Radau'use

7: scipy.integrate.solve_ivp 

from scipy.integrate import solve_ivp, odeint # Robertsonproblem（Stiff ODE） def robertson(t, y): """ Robertsonproblem（ Stiff ODE） 3 : y1 + y2 + y3 = 1 () """ y1, y2, y3 = y dy1 = -0.04 * y1 + 1e4 * y2 * y3 dy2 = 0.04 * y1 - 1e4 * y2 * y3 - 3e7 * y2**2 dy3 = 3e7 * y2**2 return [dy1, dy2, dy3] print("=" * 60) print("scipy.integrate.solve_ivp: Stiff ODE（Robertsonproblem）") print("=" * 60) # initial value y0 = [1.0, 0.0, 0.0] t_span = (0, 1e5) # 0 from 100,000 t_eval = np.logspace(-6, 5, 200) # # method1: RK45（、Stiff） print("\n1. RK45（Runge-Kutta）") import time start = time.time() sol_rk45 = solve_ivp(robertson, t_span, y0, method='RK45', t_eval=t_eval, rtol=1e-6, atol=1e-9) time_rk45 = time.time() - start print(f" calculation: {time_rk45:.3f} ") print(f" functioniterations: {sol_rk45.nfev}") print(f" : {sol_rk45.success}") # method2: BDF（、Stiff） print("\n2. BDF（、Stiff）") start = time.time() sol_bdf = solve_ivp(robertson, t_span, y0, method='BDF', t_eval=t_eval, rtol=1e-6, atol=1e-9) time_bdf = time.time() - start print(f" calculation: {time_bdf:.3f} ") print(f" functioniterations: {sol_bdf.nfev}") print(f" : {sol_bdf.success}") # method3: Radau（Runge-Kutta、Stiff） print("\n3. Radau（Runge-Kutta、Stiff）") start = time.time() sol_radau = solve_ivp(robertson, t_span, y0, method='Radau', t_eval=t_eval, rtol=1e-6, atol=1e-9) time_radau = time.time() - start print(f" calculation: {time_radau:.3f} ") print(f" functioniterations: {sol_radau.nfev}") print(f" : {sol_radau.success}") print(f"\ncalculation :") print(f" BDF vs RK45: {time_rk45 / time_bdf:.2f}timesfast") # fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5)) # : time evolution（BDF） ax1.semilogx(sol_bdf.t, sol_bdf.y[0], 'b-', linewidth=2, label='y₁') ax1.semilogx(sol_bdf.t, sol_bdf.y[1], 'r-', linewidth=2, label='y₂') ax1.semilogx(sol_bdf.t, sol_bdf.y[2], 'g-', linewidth=2, label='y₃') ax1.set_xlabel(' t [s]', fontsize=12) ax1.set_ylabel('', fontsize=12) ax1.set_title('Robertsonproblem （BDF）', fontsize=14) ax1.legend(fontsize=11) ax1.grid(True, alpha=0.3) # : (y1 + y2 + y3 = 1) conservation_bdf = sol_bdf.y[0] + sol_bdf.y[1] + sol_bdf.y[2] conservation_rk45 = sol_rk45.y[0] + sol_rk45.y[1] + sol_rk45.y[2] ax2.semilogx(sol_bdf.t, np.abs(conservation_bdf - 1), 'b-', linewidth=2, label='BDF') ax2.semilogx(sol_rk45.t, np.abs(conservation_rk45 - 1), 'r--', linewidth=2, label='RK45') ax2.set_xlabel(' t [s]', fontsize=12) ax2.set_ylabel('|y₁ + y₂ + y₃ - 1|', fontsize=12) ax2.set_title(' ', fontsize=14) ax2.legend(fontsize=11) ax2.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('stiff_ode_robertson.png', dpi=150, bbox_inches='tight') plt.show() # print(f"\n t = {sol_bdf.t[-1]:.2e} :") print(f" y₁ = {sol_bdf.y[0][-1]:.6e}") print(f" y₂ = {sol_bdf.y[1][-1]:.6e}") print(f" y₃ = {sol_bdf.y[2][-1]:.6e}") print(f" = {np.sum(sol_bdf.y[:, -1]):.10f} (: 1.0)") # : （Lotka-Volterraequation） print("\n" + "=" * 60) print("application: Lotka-Volterraequation（-）") print("=" * 60) def lotka_volterra(t, y, alpha=1.0, beta=0.1, gamma=1.5, delta=0.075): """ Lotka-Volterraequation y[0]: y[1]: """ x, y_prey = y dx = alpha * x - beta * x * y_prey dy = -gamma * y_prey + delta * x * y_prey return [dx, dy] y0 = [40, 9] # t_span = (0, 50) t_eval = np.linspace(0, 50, 500) sol_lv = solve_ivp(lotka_volterra, t_span, y0, method='RK45', t_eval=t_eval, dense_output=True) # fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5)) # ax1.plot(sol_lv.t, sol_lv.y[0], 'b-', linewidth=2, label='') ax1.plot(sol_lv.t, sol_lv.y[1], 'r-', linewidth=2, label='') ax1.set_xlabel(' t', fontsize=12) ax1.set_ylabel('', fontsize=12) ax1.set_title('Lotka-Volterraequation: ', fontsize=14) ax1.legend(fontsize=11) ax1.grid(True, alpha=0.3) # ax2.plot(sol_lv.y[0], sol_lv.y[1], 'g-', linewidth=2) ax2.plot(sol_lv.y[0][0], sol_lv.y[1][0], 'ko', markersize=10, label='initial value') ax2.set_xlabel('', fontsize=12) ax2.set_ylabel('', fontsize=12) ax2.set_title('（）', fontsize=14) ax2.legend(fontsize=11) ax2.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('lotka_volterra.png', dpi=150, bbox_inches='tight') plt.show() print("\n （ ）")
============================================================ scipy.integrate.solve_ivp: Stiff ODE（Robertsonproblem） ============================================================ 1. RK45（Runge-Kutta） calculation: 1.234 functioniterations: 8542 : True 2. BDF（、Stiff） calculation: 0.087 functioniterations: 542 : True 3. Radau（Runge-Kutta、Stiff） calculation: 0.123 functioniterations: 687 : True calculation : BDF vs RK45: 14.18timesfast t = 1.00e+05 : y₁ = 7.158272e-01 y₂ = 9.185535e-06 y₃ = 2.841636e-01 = 0.9999999999 (: 1.0) ============================================================ application: Lotka-Volterraequation（-） ============================================================ （ ）

🏋️ Exercises

1: Euler stable

Stiffequation \( dy/dt = -100y \), \( y(0) = 1 \) forward Euler method 。 \( h = 0.01, 0.02, 0.03 \) calculation、 from stable 。

2: RK4 implementation

ODERK4 、 and :

\[ \frac{dy}{dt} = t^2 + y, \quad y(0) = 1 \]

interval \([0, 1]\)、 \( h = 0.1 \) calculation、error。

3: Solution of the equation

equation \( dy/dt = r y (1 - y/K) \) 、 \( r \)（） and \( K \)（） and 。\( r = 0.5, 1.0, 2.0 \)、\( K = 1.0 \) 。

4: ODE

do ODE:

\[ \begin{cases} \frac{dx}{dt} = v \\ \frac{dv}{dt} = -2\zeta\omega_0 v - \omega_0^2 x \end{cases} \]

RK4 。: \( \omega_0 = 2\pi \)、\( \zeta = 0.1 \)（）、initial value: \( x(0) = 1, v(0) = 0 \)。 \((x, v)\) 。

Exercise 5: Application to Materials Science

equation（1、）:

\[ \frac{dT_i}{dt} = \alpha \frac{T_{i+1} - 2T_i + T_{i-1}}{\Delta x^2} \]

ODE and 。: \( T_0 = 100 \)℃、\( T_N = 0 \)℃、: 0℃。scipy.integrate.solve_ivpuse、time evolution。

Summary

、equation :

equation 、time evolution analysis 。materials science 、chemical reaction kinetics、heat conduction etc. 。 、SciPycalculation and 、suitable・・。

← Chapter 3Chapter 5 →

Disclaimer