Ordinary Differential Equations (12): Boundary Value Problems
Chen Kai BOSS
  2019-06-03 10:45:00 2019-06-03 10:45      2.9k Words  18 Mins  

Initial value problems give us complete information at the starting point, while boundary value problems provide partial information at both ends — like knowing the departure and destination points without knowing the complete path in between. These problems are extremely common in engineering: beam deflection, steady-state heat conduction, and eigenstate problems in the Schr ö dinger equation. This chapter explores in depth the theory and numerical methods for solving boundary value problems.

What Are Boundary Value Problems?

From Initial Value Problems to Boundary Value Problems

Recall the Initial Value Problem (IVP):

All conditions are given at the same point.

Boundary Value Problem (BVP):Conditions are given at different points: the left and right boundaries.

Why Are Boundary Value Problems Harder?

For initial value problems, the existence and uniqueness theorem (Picard's theorem) guarantees the existence and uniqueness of solutions. But boundary value problems may have:

  1. No solution: Incompatible conditions
  2. Infinitely many solutions: Conditions insufficient to determine a unique solution
  3. A unique solution: This is the case we hope for

Example 1: Consider

  • BVP1:,
    • Solution:, anysatisfies this, infinitely many solutions
  • BVP2:,
    • Need, contradiction, no solution
  • BVP3:,
    • Unique solution

Types of Boundary Conditions

First kind (Dirichlet): Given function values

Second kind (Neumann): Given derivative values

Third kind (Robin/Mixed): Given linear combinations of function values and derivatives

Linear Boundary Value Problems

Standard Form

Second-order linear BVP:

Solution Structure

The solution of a linear BVP can be expressed as:whereis a particular solution, and,are fundamental solutions of the homogeneous equation.

Superposition Method

Method:

  1. Solve two initial value problems:
    • Problem A:,,
    • Problem B:,,
  2. Solution is
  3. Determinefrom boundary condition:

(Requires)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp, solve_bvp
from scipy.optimize import brentq, fsolve

def superposition_method(p, q, r, a, b, alpha, beta):
    """
    Superposition method for linear BVP
    y'' + p(x)y' + q(x)y = r(x)
    y(a) = alpha, y(b) = beta
    """
    # Convert to first-order system
    def system_A(x, Y):
        y, dy = Y
        return [dy, r(x) - p(x)*dy - q(x)*y]
    
    def system_B(x, Y):
        y, dy = Y
        return [dy, -p(x)*dy - q(x)*y]
    
    # Solve problem A
    sol_A = solve_ivp(system_A, [a, b], [alpha, 0], dense_output=True)
    
    # Solve problem B
    sol_B = solve_ivp(system_B, [a, b], [0, 1], dense_output=True)
    
    # Compute constant c
    yA_b = sol_A.sol(b)[0]
    yB_b = sol_B.sol(b)[0]
    
    if abs(yB_b) < 1e-10:
        raise ValueError("Method fails: y_B(b) ≈ 0")
    
    c = (beta - yA_b) / yB_b
    
    # Construct solution
    def solution(x):
        return sol_A.sol(x)[0] + c * sol_B.sol(x)[0]
    
    return solution

# Example: y'' - y = -x, y(0) = 0, y(1) = 1
p = lambda x: 0
q = lambda x: -1
r = lambda x: -x
a, b = 0, 1
alpha, beta = 0, 1

sol = superposition_method(p, q, r, a, b, alpha, beta)

x_plot = np.linspace(a, b, 100)
y_num = np.array([sol(x) for x in x_plot])
y_exact = x_plot  # Exact solution is y = x

plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.plot(x_plot, y_num, 'b-', linewidth=2, label='Numerical')
plt.plot(x_plot, y_exact, 'r--', linewidth=2, label='Exact')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Superposition Method')
plt.legend()
plt.grid(True, alpha=0.3)

plt.subplot(122)
plt.plot(x_plot, abs(y_num - y_exact), 'k-')
plt.xlabel('x')
plt.ylabel('|Error|')
plt.title('Error')
plt.grid(True, alpha=0.3)
plt.yscale('log')

plt.tight_layout()
plt.show()

Shooting Method

Basic Idea

The Shooting Method transforms a BVP into an IVP:

  1. Guess the missing initial condition
  2. Solve the initial value problem to the right boundary
  3. Check if the right boundary condition is satisfied
  4. If not, adjust the guess and repeat

Like shooting: adjust the gun's angle (initial slope) until hitting the target (satisfying the right boundary condition).

Shooting Method for Nonlinear BVPs

Consider:We don't know.

Algorithm:

  1. Define residual function - whereis the solution of IVP with,$y'(a) = sF(s) = 0$(using Newton's method or bisection)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
def shooting_method(f, a, b, alpha, beta, s_guess, tol=1e-8, max_iter=50):
    """
    Shooting method for nonlinear BVP
    y'' = f(x, y, y')
    y(a) = alpha, y(b) = beta
    
    Returns: x array, y array
    """
    def residual(s):
        """Given initial slope s, compute boundary residual"""
        def system(x, Y):
            y, dy = Y
            return [dy, f(x, y, dy)]
        
        sol = solve_ivp(system, [a, b], [alpha, s], dense_output=True)
        return sol.sol(b)[0] - beta
    
    # Solve using Brent's method or secant method
    try:
        # Try bisection first (need to find interval containing root)
        s_low, s_high = -10, 10
        while residual(s_low) * residual(s_high) > 0:
            s_low *= 2
            s_high *= 2
        s_opt = brentq(residual, s_low, s_high, xtol=tol)
    except:
        # If bisection fails, use Newton's method
        s_opt = fsolve(residual, s_guess, full_output=False)[0]
    
    # Compute final solution with optimal s
    def system(x, Y):
        y, dy = Y
        return [dy, f(x, y, dy)]
    
    sol = solve_ivp(system, [a, b], [alpha, s_opt], dense_output=True)
    
    x_out = np.linspace(a, b, 100)
    y_out = sol.sol(x_out)[0]
    
    return x_out, y_out, s_opt

# Example: Nonlinear BVP
# y'' = -e^y, y(0) = 0, y(1) = 0
# This is a form of the Bratu problem
f = lambda x, y, dy: -np.exp(y)
a, b = 0, 1
alpha, beta = 0, 0

x, y, s = shooting_method(f, a, b, alpha, beta, s_guess=0)

plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.plot(x, y, 'b-', linewidth=2)
plt.xlabel('x')
plt.ylabel('y')
plt.title(f'Shooting Method Solution: Initial Slope s = {s:.6f}')
plt.grid(True, alpha=0.3)

# Verify boundary conditions
print(f"y(0) = {y[0]:.6f} (should be {alpha})")
print(f"y(1) = {y[-1]:.6f} (should be {beta})")

# Visualize the shooting process
s_values = np.linspace(-2, 2, 9)
plt.subplot(122)
for s_try in s_values:
    def system(x, Y):
        return [Y[1], -np.exp(Y[0])]
    sol = solve_ivp(system, [a, b], [alpha, s_try], dense_output=True)
    x_plot = np.linspace(a, b, 50)
    y_plot = sol.sol(x_plot)[0]
    plt.plot(x_plot, y_plot, '-', alpha=0.5, label=f's={s_try:.1f}')

plt.axhline(y=beta, color='r', linestyle='--', label='Target y(1)=0')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Shooting Process: Trajectories for Different Initial Slopes')
plt.legend(fontsize=8)
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('shooting_method.png', dpi=150, bbox_inches='tight')
plt.show()

Sensitivity of the Shooting Method

The shooting method can sometimes be very sensitive: the solution's dependence on initial slope can be exponential. Consider:The general solution is. At,while.

Small initial errors get amplified by about

Multiple Shooting Method

Idea: Divide the interval into multiple segments, shoot independently on each segment, then match at internal nodes.

This reduces exponential amplification of errors.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
def multiple_shooting(f, a, b, alpha, beta, n_segments=4):
    """
    Multiple shooting method
    Divides [a,b] into n_segments segments
    """
    # Nodes
    nodes = np.linspace(a, b, n_segments + 1)
    h = nodes[1] - nodes[0]
    
    # Initial guess: linear interpolation
    y_guess = np.linspace(alpha, beta, n_segments + 1)
    s_guess = np.ones(n_segments) * (beta - alpha) / (b - a)
    
    def residuals(params):
        """
        params = [y_1, ..., y_{n-1}, s_0, s_1, ..., s_{n-1}]
        y_0 = alpha, y_n = beta are fixed
        """
        y_interior = params[:n_segments-1]
        s_all = params[n_segments-1:]
        
        y_nodes = np.concatenate([[alpha], y_interior, [beta]])
        
        res = []
        
        for i in range(n_segments):
            x_left, x_right = nodes[i], nodes[i+1]
            y_left, s_left = y_nodes[i], s_all[i]
            
            def system(x, Y):
                return [Y[1], f(x, Y[0], Y[1])]
            
            sol = solve_ivp(system, [x_left, x_right], [y_left, s_left])
            y_right_computed = sol.y[0, -1]
            s_right_computed = sol.y[1, -1]
            
            # Continuity conditions
            res.append(y_right_computed - y_nodes[i+1])
            if i < n_segments - 1:
                res.append(s_right_computed - s_all[i+1])  # Slope continuity
        
        return res
    
    # Initial parameters
    initial_params = np.concatenate([y_guess[1:-1], s_guess])
    
    # Solve
    solution = fsolve(residuals, initial_params)
    
    # Reconstruct solution
    y_interior = solution[:n_segments-1]
    s_all = solution[n_segments-1:]
    y_nodes = np.concatenate([[alpha], y_interior, [beta]])
    
    # Solve on each segment and merge
    x_full, y_full = [], []
    for i in range(n_segments):
        x_left, x_right = nodes[i], nodes[i+1]
        
        def system(x, Y):
            return [Y[1], f(x, Y[0], Y[1])]
        
        sol = solve_ivp(system, [x_left, x_right], [y_nodes[i], s_all[i]],
                       dense_output=True)
        
        x_seg = np.linspace(x_left, x_right, 25)
        y_seg = sol.sol(x_seg)[0]
        
        x_full.extend(x_seg[:-1] if i < n_segments-1 else x_seg)
        y_full.extend(y_seg[:-1] if i < n_segments-1 else y_seg)
    
    return np.array(x_full), np.array(y_full)

Finite Difference Method

Basic Idea

The Finite Difference Method directly discretizes the differential equation, transforming the BVP into a system of algebraic equations.

Discretizing derivatives:

Finite Differences for Linear BVPs

Consider: Discretize at nodes(,):Rearranging:

Matrix Form

This is a tridiagonal linear system:where: - - - -(needs boundary term corrections)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
def finite_difference_linear(p, q, r, a, b, alpha, beta, N):
    """
    Finite difference method for linear BVP
    y'' + p(x)y' + q(x)y = r(x)
    y(a) = alpha, y(b) = beta
    """
    h = (b - a) / N
    x = np.linspace(a, b, N + 1)
    
    # Interior nodes
    x_interior = x[1:-1]
    n = len(x_interior)
    
    # Construct tridiagonal matrix
    A = np.zeros((n, n))
    b_vec = np.zeros(n)
    
    for i in range(n):
        xi = x_interior[i]
        pi, qi, ri = p(xi), q(xi), r(xi)
        
        # Diagonal
        A[i, i] = -2 + h**2 * qi
        
        # Lower diagonal
        if i > 0:
            A[i, i-1] = 1 - h * pi / 2
        
        # Upper diagonal
        if i < n - 1:
            A[i, i+1] = 1 + h * pi / 2
        
        # Right-hand side
        b_vec[i] = h**2 * ri
        
        # Boundary corrections
        if i == 0:
            b_vec[i] -= (1 - h * pi / 2) * alpha
        if i == n - 1:
            b_vec[i] -= (1 + h * pi / 2) * beta
    
    # Solve linear system
    y_interior = np.linalg.solve(A, b_vec)
    
    # Add boundary values
    y = np.concatenate([[alpha], y_interior, [beta]])
    
    return x, y

# Example: y'' - y = -x, y(0) = 0, y(1) = 1
p = lambda x: 0
q = lambda x: -1
r = lambda x: -x
a, b = 0, 1
alpha, beta = 0, 1

# Different grid densities
N_values = [10, 20, 40, 80]

plt.figure(figsize=(14, 5))

plt.subplot(121)
for N in N_values:
    x, y = finite_difference_linear(p, q, r, a, b, alpha, beta, N)
    plt.plot(x, y, 'o-', markersize=3, label=f'N = {N}')

# Exact solution
x_exact = np.linspace(a, b, 200)
y_exact = x_exact
plt.plot(x_exact, y_exact, 'k-', linewidth=2, label='Exact')

plt.xlabel('x')
plt.ylabel('y')
plt.title('Finite Difference Method')
plt.legend()
plt.grid(True, alpha=0.3)

# Convergence analysis
plt.subplot(122)
errors = []
for N in [10, 20, 40, 80, 160, 320]:
    x, y = finite_difference_linear(p, q, r, a, b, alpha, beta, N)
    y_exact_at_nodes = x
    error = np.max(np.abs(y - y_exact_at_nodes))
    errors.append(error)

h_values = (b - a) / np.array([10, 20, 40, 80, 160, 320])
plt.loglog(h_values, errors, 'bo-', label='Maximum error')
plt.loglog(h_values, h_values**2, 'r--', label='O(h ²) reference')
plt.xlabel('Step size h')
plt.ylabel('Maximum error')
plt.title('Convergence Order Verification')
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('finite_difference.png', dpi=150, bbox_inches='tight')
plt.show()

Finite Differences for Nonlinear BVPs

For nonlinear equations, discretization yields a nonlinear algebraic system that requires Newton's method to solve.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
def finite_difference_nonlinear(f, df_dy, df_ddy, a, b, alpha, beta, N, 
                                 tol=1e-10, max_iter=100):
    """
    Finite difference method for nonlinear BVP (Newton iteration)
    y'' = f(x, y, y')
    y(a) = alpha, y(b) = beta
    """
    h = (b - a) / N
    x = np.linspace(a, b, N + 1)
    n = N - 1  # Number of interior nodes
    
    # Initial guess: linear interpolation
    y = np.linspace(alpha, beta, N + 1)
    
    for iteration in range(max_iter):
        # Construct residual vector
        F = np.zeros(n)
        
        for i in range(1, N):
            xi = x[i]
            yi = y[i]
            dyi = (y[i+1] - y[i-1]) / (2*h)
            d2yi = (y[i+1] - 2*y[i] + y[i-1]) / h**2
            
            F[i-1] = d2yi - f(xi, yi, dyi)
        
        # Check convergence
        if np.max(np.abs(F)) < tol:
            break
        
        # Construct Jacobian matrix
        J = np.zeros((n, n))
        
        for i in range(1, N):
            xi = x[i]
            yi = y[i]
            dyi = (y[i+1] - y[i-1]) / (2*h)
            
            # Partial derivative with respect to y[i-1]
            if i > 1:
                J[i-1, i-2] = 1/h**2 + df_ddy(xi, yi, dyi) / (2*h)
            
            # Partial derivative with respect to y[i]
            J[i-1, i-1] = -2/h**2 - df_dy(xi, yi, dyi)
            
            # Partial derivative with respect to y[i+1]
            if i < N - 1:
                J[i-1, i] = 1/h**2 - df_ddy(xi, yi, dyi) / (2*h)
        
        # Newton update
        dy = np.linalg.solve(J, -F)
        y[1:-1] += dy
    
    return x, y

Eigenvalue Problems

Sturm-Liouville Problems

Standard form: Hereis the eigenvalue andis the eigenfunction.

Finding Eigenvalues with Finite Differences

After discretization, we get a generalized eigenvalue problem.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
def sturm_liouville_eigenvalues(p, q, w, a, b, N, num_eig=5):
    """
    Sturm-Liouville eigenvalue problem
    -(p(x)y')' + q(x)y = λ w(x)y
    y(a) = y(b) = 0
    """
    h = (b - a) / N
    x = np.linspace(a, b, N + 1)
    n = N - 1
    
    # Construct matrix A (from -(py')' + qy)
    A = np.zeros((n, n))
    B = np.zeros((n, n))  # Mass matrix (from wy)
    
    for i in range(n):
        xi = x[i + 1]
        
        # Main diagonal
        p_plus = p(xi + h/2)
        p_minus = p(xi - h/2)
        A[i, i] = (p_plus + p_minus) / h**2 + q(xi)
        B[i, i] = w(xi)
        
        # Off-diagonals
        if i > 0:
            A[i, i-1] = -p_minus / h**2
        if i < n - 1:
            A[i, i+1] = -p_plus / h**2
    
    # Solve generalized eigenvalue problem
    from scipy.linalg import eigh
    eigenvalues, eigenvectors = eigh(A, B)
    
    # Reconstruct eigenfunctions (add boundary value 0)
    eigenfunctions = []
    for j in range(num_eig):
        ef = np.concatenate([[0], eigenvectors[:, j], [0]])
        # Normalize
        ef = ef / np.sqrt(np.sum(ef**2 * h))
        eigenfunctions.append(ef)
    
    return x, eigenvalues[:num_eig], eigenfunctions

# Example: Simple Sturm-Liouville problem
# -y'' = λ y, y(0) = y(π) = 0
# Exact eigenvalues: λ_n = n ², eigenfunctions: sin(nx)
p = lambda x: 1
q = lambda x: 0
w = lambda x: 1
a, b = 0, np.pi
N = 100

x, eigenvalues, eigenfunctions = sturm_liouville_eigenvalues(p, q, w, a, b, N, num_eig=6)

# Compare with exact values
exact_eigenvalues = np.array([1, 4, 9, 16, 25, 36])

print("Sturm-Liouville Eigenvalues:")
print("Numerical\tExact\t\tRelative Error")
for i, (num, exact) in enumerate(zip(eigenvalues, exact_eigenvalues)):
    rel_err = abs(num - exact) / exact
    print(f"{num:.6f}\t{exact}\t\t{rel_err:.2e}")

Practical Applications

Beam Deflection

The deflectionof a simply supported beam satisfies:whereis Young's modulus,is the moment of inertia, andis the load distribution.

Boundary conditions (simply supported):

Steady-State Heat Conduction

One-dimensional steady-state heat conduction:Boundary conditions can be: - Temperature specified (Dirichlet) - Heat flux specified (Neumann) - Convective boundary (Robin)

Quantum Mechanics: Schr ö dinger Equation

One-dimensional time-independent Schr ö dinger equation:This is an eigenvalue problem!

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
def schrodinger_1d(V, x_min, x_max, N, num_states=5, hbar=1, m=1):
    """
    One-dimensional time-independent Schr ö dinger equation
    Returns energy levels and wavefunctions
    """
    h = (x_max - x_min) / N
    x = np.linspace(x_min, x_max, N + 1)
    n = N - 1  # Interior points
    
    # Hamiltonian matrix
    H = np.zeros((n, n))
    
    kinetic_coeff = hbar**2 / (2 * m * h**2)
    
    for i in range(n):
        xi = x[i + 1]
        
        # Kinetic energy term (-ℏ²/(2m) d ²/dx ²)
        H[i, i] = 2 * kinetic_coeff + V(xi)
        if i > 0:
            H[i, i-1] = -kinetic_coeff
        if i < n - 1:
            H[i, i+1] = -kinetic_coeff
    
    # Solve eigenvalue problem
    energies, wavefunctions = np.linalg.eigh(H)
    
    # Reconstruct wavefunctions (add boundary value 0)
    psi = []
    for j in range(num_states):
        wf = np.concatenate([[0], wavefunctions[:, j], [0]])
        # Normalize
        wf = wf / np.sqrt(np.sum(wf**2) * h)
        psi.append(wf)
    
    return x, energies[:num_states], psi

# Example 1: Infinite square well
# V = 0 inside, V = ∞ at boundaries
# Exact energy levels: E_n = n ²π²ℏ²/(2mL ²)
V_well = lambda x: 0
L = 1
x, E_well, psi_well = schrodinger_1d(V_well, 0, L, N=200, num_states=4)

# Example 2: Harmonic oscillator potential
# V = (1/2)kx ² = (1/2)m ω² x ²
omega = 1
V_harmonic = lambda x: 0.5 * omega**2 * x**2
x_harm, E_harm, psi_harm = schrodinger_1d(V_harmonic, -5, 5, N=200, num_states=4)
# Exact energy levels: E_n = (n + 1/2)ℏω

Using SciPy to Solve BVPs

SciPy provides the solve_bvp function:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
from scipy.integrate import solve_bvp

# Example: y'' = -e^y, y(0) = 0, y(1) = 0
def fun(x, y):
    """Convert second-order equation to first-order system"""
    return np.vstack((y[1], -np.exp(y[0])))

def bc(ya, yb):
    """Boundary conditions: ya is value at x=a, yb is value at x=b"""
    return np.array([ya[0] - 0,   # y(0) = 0
                     yb[0] - 0])  # y(1) = 0

# Initial mesh and guess
x = np.linspace(0, 1, 10)
y = np.zeros((2, x.size))  # Initial guess

# Solve
sol = solve_bvp(fun, bc, x, y)

if sol.success:
    x_plot = np.linspace(0, 1, 100)
    y_plot = sol.sol(x_plot)[0]
    
    plt.figure(figsize=(10, 5))
    plt.plot(x_plot, y_plot, 'b-', linewidth=2)
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('solve_bvp Solution of Bratu Problem')
    plt.grid(True, alpha=0.3)
    plt.show()
else:
    print("Solver failed:", sol.message)

Summary

In this chapter, we learned various numerical methods for boundary value problems:

Method Advantages Disadvantages Suitable For
Superposition Simple, uses existing IVP solvers Linear problems only Linear BVPs
Shooting Flexible, handles nonlinear Can be unstable General BVPs
Multiple Shooting More stable Complex implementation Sensitive BVPs
Finite Difference Direct, easy to understand Requires linear system solve Regular domains
Finite Element Handles complex geometry Theoretically complex Engineering problems

Practical advice:

  1. First try scipy.integrate.solve_bvp
  2. If it fails, check your initial guess
  3. For sensitive problems, consider multiple shooting or finite differences
  4. For eigenvalue problems, finite differences are usually sufficient

In the next chapter, we will enter the world of partial differential equations— what happens when the unknown function depends on multiple variables?

  • Post title：Ordinary Differential Equations (12): Boundary Value Problems
  • Post author：Chen Kai
  • Create time：2019-06-03 10:45:00
  • Post link：https://www.chenk.top/ode-chapter-12-boundary-value-problems/
  • Copyright Notice：All articles in this blog are licensed under BY-NC-SA unless stating additionally.
 Comments