<h1 style="color:purple;font-size:3em"> Boundary Value Problem </h1>

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:

<h2 style="color:purple;font-size:2em"> Shooting Method </h2>

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$ and for 5, $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 
$$


```python
# importing necessary modules
import numpy as np
import matplotlib.pyplot as plt
import scipy
```


```python
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
```

Lets check the result with two initial guess of $\alpha$, $\alpha_1 = 1$ and $\alpha_2 = -10$


```python
desolve(1)[-1]    # last term of solution with alpha = 1
```




    1.6346151809563805




```python
desolve(-10)[-1]   # last term of solution with alpha = 5
```




    -1.0650168118209011



This makes a good guess, since both have opposite sign.

Now, using bisection method to find the value of alpha : 


```python
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 =  -4.5
    alpha =  -7.25
    alpha =  -5.875
    alpha =  -5.1875
    alpha =  -5.53125
    alpha =  -5.703125
    alpha =  -5.6171875
    alpha =  -5.66015625
    alpha =  -5.681640625
    alpha =  -5.6708984375
    alpha =  -5.66552734375
    alpha =  -5.662841796875
    alpha =  -5.6614990234375
    alpha =  -5.66082763671875
    alpha =  -5.660491943359375
    alpha =  -5.6603240966796875
    alpha =  -5.660408020019531
    alpha =  -5.660449981689453
    alpha =  -5.660470962524414
    alpha =  -5.660460472106934
    alpha =  -5.660455226898193
    alpha =  -5.660452604293823
    alpha =  -5.660451292991638
    alpha =  -5.660450637340546
    alpha =  -5.660450965166092



```python
alpha
```




    -5.660451129078865




```python
#  Now we  have alpha, we can find the solution
y = desolve(alpha)
y[-1]
```




    3.284482948638645e-09




```python
# 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]
```




    (1.0, -1.3877787807814457e-17)




```python
# 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()
```

### Newton Raphson Shooting Method

Remember 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
$$


```python
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
```

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
$$



```python
# 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 =  1



```python
alpha
```




    -5.660451142461925




```python
# Now, we can calculate value of y for this alpha
ynr = desolve(alpha)
ynr[-1]
```




    -1.951563910473908e-15




```python
# 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()
```

<h2 style="color:purple;font-size:2em"> Linear Shooting Method </h2>

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


```python
# 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
```


```python
# 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
```


```python
# 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
```


```python
# storing arrays 
u = desolve1()
v = desolve2()
```


```python
# checking if v(x_n) = 0:
v[-1]
```




    0.24542109025247982



Now the solution is given by
$$
y(x) = u(x) + \frac{b - u(x_n)}{v(x_n)}v(x)
$$


```python
ylin = u + (-u[-1]/v[-1])*v
ylin[-1]
```




    0.0




```python
plt.figure(figsize=(10, 8))
plt.plot(x, ylin, lw=5, label="linear shooting")
plt.plot(x, ytrue, "--", lw=5, label="exact")
plt.legend()
```


```python

```


```python

```

9c60b71e441410fc47885e511af0d7b9d22f0594
ipynb
Jupyter Notebook
differential_equations/boundary_value_problem.ipynb
phy6boy/algorithms
affe3abcf3cd5261aa283f2ae5d1c1dc3165f8de
MIT
differential_equations/boundary_value_problem.ipynb
phy6boy/algorithms
affe3abcf3cd5261aa283f2ae5d1c1dc3165f8de
MIT
differential_equations/boundary_value_problem.ipynb
phy6boy/algorithms
affe3abcf3cd5261aa283f2ae5d1c1dc3165f8de
MIT
Qwen/Qwen-72B
1. YES
2. YES