<jupyter_start><jupyter_text>Boundary Value Problem  As opposed to initial value problems, [Boundary value problems](https://en.wikipedia.org/wiki/Boundary_value_problem) is a class of differential equation together with a set of additional constraints (boundary conditions). A BVP with second order D.E has the form : $$\frac{d^2y}{dx^2}\ =\ f\left(x, y(x), y^{\prime}(x)\right),\quad y(x_0) = a,\quad y(x_n) = b$$Unlike initial value problems, a BVP has no unique solution. Some of theqniques used for finding solution are described here:  Shooting Method In this approach, we convert the BVP to Initial Value Problem(IVP). So, if we take $u_\alpha$ to be the solution of the IVP, equivalent to the previous BVP, then it has the form :$$\frac{d^2u_\alpha}{dx^2}\ =\ f\left(x, u_\alpha(x), u_\alpha^{\prime}(x)\right),\quad u_\alpha(x_0) = a,\quad u^\prime_\alpha(x_0) = \alpha$$Here, $\alpha$ is choosen with an educated guess. Because of this guess, the value of this function $u_\alpha(x)$ at $x=x_n$ will not be same as $y(x_n)$. Let us take the difference as $g(\alpha)$, since the this depends on $\alpha$. Thus we have : $$g(\alpha)\ =\ u_\alpha(x_n) - y(x_n)\ =\ u_\alpha(x_n) - b \qquad\qquad \because y(x_n) = b$$So, inorder to solve this set of equations, one has to find zero of equation above. Here are some expamples using different methods:  Bisection Shooting Method**Question :** solve the D.E :$$y^{\prime\prime} + 4y^\prime\ =\ 5x, \quad y(0)=1, \quad y(1)=0$$**Solution :**Lets convert this to a initial value problem. $$u_\alpha^{\prime\prime} + 4u_\alpha^\prime\ =\ 5x, \quad u_\alpha(0)=1, \quad u^\prime_\alpha(0)=\alpha$$lets pick initial values of $\alpha$ are 0 and 5, so that when we take it as 0 and solve, $u_\alpha(1)0$. Lets try implimenting a runge kutta method for finding D.E  Runge Kutta the equation :$$y^{\prime\prime} + 4y^\prime\ =\ 5x$$can be split into two differential equations by substituting for $y^\prime = u$, so that the equation becomes\begin{align}\frac{dy}{dx} &= u \\[0.5em]\frac{du}{dx} &= 5x-4u\end{align}with constraints$$y(0) = 1, \quad u(0) = \alpha $$<jupyter_code># importing necessary modules
import numpy as np
import matplotlib.pyplot as plt
import scipy
N = 200                      # Number of samples
x = np.linspace(0, 1, N)     # Creating the array x
dx = x[1] - x[0]             # Step size


# defining f(x, u) for derivatve du/dx = f(x, u)
def f(x, u):
    return 5*x - 4*u

# The following function takes a value of alpha and find the solution of the equation
# y'' + 4y' = 5x using Runge-kutta method

def desolve(alpha):
    y = np.zeros(N, dtype=float) # Creating array to store y
    u = np.zeros(N, dtype=float) # Creating array to store u
    
    # setting initial values
    y[0] = 1
    u[0] = alpha
    
    # Starting R-K loop
    for i in range(1, N):
        # k_x represents for y(x) and l_x for u(x)
        k1 = u[i-1]
        l1 = f(x[i-1], u[i-1])
        k2 = u[i-1] + (1/2)*l1*dx
        l2 = f(x[i-1] + (1/2)*dx, u[i-1] + (1/2)*l1*dx)
        k3 = u[i-1] + (1/2)*l2*dx
        l3 = f(x[i-1] + (1/2)*dx, u[i-1] + (1/2)*l2*dx)
        k4 = u[i-1] + l3*dx
        l4 = f(x[i-1] + dx, u[i-1] + l3*dx)
        
        # calculating the next y and u values
        y[i] = y[i-1] + (1/6)*dx*(k1 + 2*k2 + 2*k3 + k4)
        u[i] = u[i-1] + (1/6)*dx*(l1 + 2*l2 + 2*l3 + l4)
    # Now, the array y stores the result
    return y<jupyter_output><empty_output><jupyter_text>Lets check the result with two initial guess of $\alpha$, $\alpha_1 = 1$ and $\alpha_2 = -10$<jupyter_code>desolve(1)[-1]    # last term of solution with alpha = 1
desolve(-10)[-1]   # last term of solution with alpha = 5<jupyter_output><empty_output><jupyter_text>This makes a good guess, since both have opposite sign.Now, using bisection method to find the value of alpha :<jupyter_code>a = 1
b = -10
# bisection
alpha = (a + b)/2
# final value of solution at alpha
fa = desolve(alpha)[-1]
while abs(fa)>1e-8:
    print("alpha = ", alpha)
    # If the last value of y is positive, then update a, else updtae b
    if fa > 0 :
        a = alpha
    else:
        b = alpha
    # update alpha and fa
    alpha = (a + b)/2
    fa = desolve(alpha)[-1]

alpha
#  Now we  have alpha, we can find the solution
y = desolve(alpha)
y[-1]
# True solution find using a Computer Algebra System, for comparison
N = 200
x = np.linspace(0, 1, N)
ytrue = (5/8)*x**2 - (5/16)*x - (1/64)*(25*np.exp(4) + 59)/(np.exp(4) - 1) +\
        (21/16)*np.exp(-4*x + 4)/(np.exp(4) - 1) + (5/64)
ytrue[0], ytrue[-1]
# visualizing the solutions
plt.figure(figsize=(10, 8))
plt.plot(x, y,"-", lw=5, label="calculated")
plt.plot(x, ytrue,"--",lw=5, label="exact")
plt.plot(x, desolve(1), "--", label="alpha=1")
plt.plot(x, desolve(-10), "--", label="alpha=-10")
plt.legend()<jupyter_output><empty_output><jupyter_text>Newton Raphson Shooting MethodRemember we have to solve IVP: $$\frac{d^2u_\alpha}{dx^2}\ =\ f\left(x, u_\alpha(x), u_\alpha^{\prime}(x)\right),\quad u_\alpha(x_0) = a,\quad u^\prime_\alpha(x_0) = \alpha$$with $g(\alpha)$ is given by:$$g(\alpha)\ =\ u_\alpha(x_n) - y(x_n)\ =\ u_\alpha(x_n) - b \qquad\qquad \because y(x_n) = b$$Here, to find root of equation $g(\alpha) = 0$, we use Newton-Raphson method:$$\alpha_{i+1} = \alpha_i - \frac{g(\alpha_i)}{g^\prime(\alpha_i)}$$But the problem here is, we don't have $g^\prime(\alpha)$. To calculate thats, we take derivative w.r.t $\alpha$ of the first equation:\begin{align}\frac{d}{d\alpha}\frac{d^2u_\alpha}{dx^2}\ &=\ \frac{d}{d\alpha}f\left(x, u_\alpha(x), u_\alpha^{\prime}(x)\right)\\[0.5em]\implies \frac{d^2}{dx^2}\frac{du_\alpha}{d\alpha}\ &=\ \frac{\partial f}{\partial u_\alpha}\frac{du_\alpha}{d\alpha} + \frac{\partial f}{\partial u^\prime_\alpha}\frac{du^\prime_\alpha}{d\alpha}\end{align}putting $\frac{du_\alpha}{d\alpha}\ =\ v_\alpha(x)$, we have,$$v^{\prime\prime}_\alpha\ =\ \frac{\partial f}{\partial u_\alpha}v_\alpha + \frac{\partial f}{\partial u^\prime_\alpha}v^\prime_\alpha$$initial conditions of $v_\alpha(x)$ are:$$v_\alpha(x_0) = \frac{du_\alpha(x_0)}{d\alpha} = \frac{da}{d\alpha} = 0, \qquad v^\prime_\alpha(x_0) = \frac{du^\prime_\alpha(x_0)}{d\alpha} = \frac{d\alpha}{d\alpha} = 1$$Now we have another D.E to solve, as described above. As a result, the newton raphson ormula becomes:$$\alpha_{i+1} = \alpha_i - \frac{u_{\alpha_i}(x_n)-b}{v_{\alpha_i}(x_n)}$$ So, inorder to solve the equation using newton raphson method, we have two IVP:$$\frac{d^2u_\alpha}{dx^2}\ =\ f\left(x, u_\alpha(x), u_\alpha^{\prime}(x)\right),\quad u_\alpha(x_0) = a,\quad u^\prime_\alpha(x_0) = \alpha$$$$v^{\prime\prime}_\alpha\ =\ \frac{\partial f}{\partial u_\alpha}v_\alpha + \frac{\partial f}{\partial u^\prime_\alpha}v^\prime_\alpha, \quad v_\alpha(x_0) = 0, \quad v^\prime_\alpha(x_0) = 1$$and use this formula to better estimate $\alpha$ and repeat until convergence is achieved$$\alpha_{i+1} = \alpha_i - \frac{u_{\alpha_i}(x_n)-b}{v_{\alpha_i}(x_n)}$$ **Question :** solve the D.E :$$y^{\prime\prime} + 4y^\prime\ =\ 5x, \quad y(0)=1, \quad y(1)=0$$**soluton :** This is same as the one we encounter in bisection method, so we can use the same desolve(alpha) to find solution for first IVP. for the second case, we will create desolve2(alpha)we have The IVP:$$u_\alpha^{\prime\prime} + 4u_\alpha^\prime\ =\ 5x, \quad u_\alpha(0)=1, \quad u^\prime_\alpha(0)=\alpha$$$$\therefore f = 5x - 4u^\prime_\alpha$$Hence$$\frac{\partial f}{\partial u_\alpha} = 0, \quad \frac{\partial f}{\partial u^\prime_\alpha} = -4$$so the 2nd IVP is:$$v^{\prime\prime}_\alpha\ = -4v^\prime_\alpha, \quad v_\alpha(x_0) = 0, \quad v^\prime_\alpha(x_0) = 1$$ So, the DE to solve is, $$v^{\prime\prime} = -4v^\prime$$we make substitutio u = $v^\prime$, so that we got two coupled 1st order DE:$$\frac{dv}{dx} = u(x), \quad \texttt{and} \quad \frac{du}{dx} = -4u(x)$$with$$v(x_0) = 0 \quad \texttt{and} \quad u(x_0) = 1$$<jupyter_code>N = 200                      # Number of samples
x = np.linspace(0, 1, N)     # Creating the array x
dx = x[1] - x[0]             # Step size


# defining f2(x, u) for derivatve du/dx = f2(x, u)
def f2(x, u):
    return -4*u

# The following function takes a value of alpha and find the solution of the equation
# y'' + 4y' = 5x using Runge-kutta method

def desolve2(alpha):
    v = np.zeros(N, dtype=float) # Creating array to store y
    u = np.zeros(N, dtype=float) # Creating array to store u
    
    # setting initial values
    v[0] = 0
    u[0] = 1
    
    # Starting R-K loop
    for i in range(1, N):
        # k_x represents for v(x) and l_x for u(x)
        k1 = u[i-1]
        l1 = f2(x[i-1], u[i-1])
        k2 = u[i-1] + (1/2)*l1*dx
        l2 = f2(x[i-1] + (1/2)*dx, u[i-1] + (1/2)*l1*dx)
        k3 = u[i-1] + (1/2)*l2*dx
        l3 = f2(x[i-1] + (1/2)*dx, u[i-1] + (1/2)*l2*dx)
        k4 = u[i-1] + l3*dx
        l4 = f2(x[i-1] + dx, u[i-1] + l3*dx)
        
        # calculating the next y and u values
        v[i] = v[i-1] + (1/6)*dx*(k1 + 2*k2 + 2*k3 + k4)
        u[i] = u[i-1] + (1/6)*dx*(l1 + 2*l2 + 2*l3 + l4)
    # Now, the array y stores the result
    return v<jupyter_output><empty_output><jupyter_text>Now we have Newton-Raphson formula:$$\alpha_{i+1} = \alpha_i - \frac{u_{\alpha_i}(x_n)-b}{v_{\alpha_i}(x_n)} = \alpha_i - \frac{u_{\alpha_i}(x_n)}{v_{\alpha_i}(x_n)} \quad \because b=0$$<jupyter_code># starting guess for alpha
alpha = 1
# taking end point
fa =  desolve(alpha)[-1]
while abs(fa)>1e-8:
    print("alpha = ", alpha)
    alpha = alpha - desolve(alpha)[-1]/desolve2(alpha)[-1]
    fa = desolve(alpha)[-1]
alpha
# Now, we can calculate value of y for this alpha
ynr = desolve(alpha)
ynr[-1]
# visualizing the solutions
plt.figure(figsize=(10, 8))
plt.plot(x, ynr,"-", lw=5, label="newton-raphson")
plt.plot(x, ytrue,".",lw=5, label="exact")

plt.legend()<jupyter_output><empty_output><jupyter_text>Linear Shooting Method A linear equation has the form : $$y^{\prime\prime}(x) + P(x)y^\prime(x) + Q(x)y(x) = R(x)$$If $y_1(x)$ and $y_2(x)$ are solutions of this equation, then so is $k_1y_1(x) + k_2y_2(x)$. This linearity makes it possible to not specifiying $\alpha$ and use root finding as did in earlier two methods. cosider boundary value problem$$y^{\prime\prime} + P(x)y^\prime(x) + Q(x)y(x) = R(x), \quad y(x_0) = a, \quad y(x_n) = b$$we can construct two IVP:$$u^{\prime\prime} + P(x)u^\prime(x) + Q(x)u(x) = R(x), \quad u(x_0) = a, \quad u^\prime(x_0) = 0$$$$v^{\prime\prime} + P(x)v^\prime(x) + Q(x)v(x) = 0, \quad v(x_0) = 0, \quad v^\prime(x_0) = s$$Here, s is a parameter. Now, if we take a linear combination of these two equations,$$\omega(x) = u(x) + \lambda v(x)$$This $\omega(x)$ satisfy left boundary condition:$$\omega(x_0) = u(x_0) + \lambda v(x_0) = a$$inorder to satisfy right boundary condition, we have $$\omega(x_n) = u(x_n) + \lambda v(x_n) = b$$this gives, $$\lambda = \frac{b - u(x_n)}{v(x_n)}$$and hence the solution, is$$y(x) = u(x) + \frac{b - u(x_n)}{v(x_n)}v(x)$$ So, as a summary, we need to setup two initial value problems and use the last equation to find the solution. This does not include any root finding steps, but if $v(x_n) = 0$, then the solution to BVP may not exist. **Question :** solve the D.E :$$y^{\prime\prime} + 4y^\prime\ =\ 5x, \quad y(0)=1, \quad y(1)=0$$ **solution :** This is a linear equation. so, we set up two IVPs (with s=1):$$u^{\prime\prime} + 4u^\prime = 5x, \quad u(0) = 1, \quad u^\prime(0) = 0$$$$v^{\prime\prime} + 4v^\prime = 0, \quad v(0) = 0, \quad v^\prime(0) = 1$$ solving both is same as we done above. so we just copy and paste it here, make simple changes<jupyter_code># store x array
N = 200                      # Number of samples
x = np.linspace(0, 1, N)     # Creating the array x
dx = x[1] - x[0]             # Step size
# Non-Homogeneous part

# defining f1(x, u) for derivatve du/dx = f1(x, u)
def f1(x, u):
    return 5*x - 4*u

# The following function takes a value of alpha and find the solution of the equation
# y'' + 4y' = 5x using Runge-kutta method

def desolve1():
    y = np.zeros(N, dtype=float) # Creating array to store y
    u = np.zeros(N, dtype=float) # Creating array to store u
    
    # setting initial values
    y[0] = 1
    u[0] = 0
    
    # Starting R-K loop
    for i in range(1, N):
        # k_x represents for y(x) and l_x for u(x)
        k1 = u[i-1]
        l1 = f1(x[i-1], u[i-1])
        k2 = u[i-1] + (1/2)*l1*dx
        l2 = f1(x[i-1] + (1/2)*dx, u[i-1] + (1/2)*l1*dx)
        k3 = u[i-1] + (1/2)*l2*dx
        l3 = f1(x[i-1] + (1/2)*dx, u[i-1] + (1/2)*l2*dx)
        k4 = u[i-1] + l3*dx
        l4 = f1(x[i-1] + dx, u[i-1] + l3*dx)
        
        # calculating the next y and u values
        y[i] = y[i-1] + (1/6)*dx*(k1 + 2*k2 + 2*k3 + k4)
        u[i] = u[i-1] + (1/6)*dx*(l1 + 2*l2 + 2*l3 + l4)
    # Now, the array y stores the result
    return y
# Homogeneous part
# defining f2(x, u) for derivatve du/dx = f2(x, u)
def f2(x, u):
    return - 4*u

# The following function takes a value of alpha and find the solution of the equation
# y'' + 4y' = 5x using Runge-kutta method

def desolve2():
    y = np.zeros(N, dtype=float) # Creating array to store y
    u = np.zeros(N, dtype=float) # Creating array to store u
    
    # setting initial values
    y[0] = 0
    u[0] = 1
    
    # Starting R-K loop
    for i in range(1, N):
        # k_x represents for y(x) and l_x for u(x)
        k1 = u[i-1]
        l1 = f2(x[i-1], u[i-1])
        k2 = u[i-1] + (1/2)*l1*dx
        l2 = f2(x[i-1] + (1/2)*dx, u[i-1] + (1/2)*l1*dx)
        k3 = u[i-1] + (1/2)*l2*dx
        l3 = f2(x[i-1] + (1/2)*dx, u[i-1] + (1/2)*l2*dx)
        k4 = u[i-1] + l3*dx
        l4 = f2(x[i-1] + dx, u[i-1] + l3*dx)
        
        # calculating the next y and u values
        y[i] = y[i-1] + (1/6)*dx*(k1 + 2*k2 + 2*k3 + k4)
        u[i] = u[i-1] + (1/6)*dx*(l1 + 2*l2 + 2*l3 + l4)
    # Now, the array y stores the result
    return y
# storing arrays 
u = desolve1()
v = desolve2()
# checking if v(x_n) = 0:
v[-1]<jupyter_output><empty_output><jupyter_text>Now the solution is given by$$y(x) = u(x) + \frac{b - u(x_n)}{v(x_n)}v(x)$$<jupyter_code>ylin = u + (-u[-1]/v[-1])*v
ylin[-1]
plt.figure(figsize=(10, 8))
plt.plot(x, ylin, lw=5, label="linear shooting")
plt.plot(x, ytrue, "--", lw=5, label="exact")
plt.legend()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19124729531653523, 'fraction_numerical': 0.036574300272213305, 'mean_word_length': 3.1602787456445993, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 32, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '434750', 'license': 'MIT', 'n_tokens_mistral': 5953, 'n_tokens_neox': 5734, 'n_words': 1618, 'path': 'differential_equations/boundary_value_problem.ipynb', 'repo_name': 'phy6boy/algorithms', 'size': '14327'}