# Shooting Method- Boundary value problem starting from -1 to 1

The equation is $$\rho \frac{d \bar{u}}{dy} = -\frac{d\bar{p}}{dx} + \mu\frac{d^2\bar{u}}{dy^2}$$

with boundary condition

$$u(-1)=0$$ and $$u(1)=1$$

I am to solve it using fifth order runge-kutta felhberg approach with Matlab. To do that, I have to change the equation to initial value problem using shooting method.

So, my problems are:

1. how can I do that when the boundary is from -1
2. Please I need the matlab code to do it. Thanks

In order to use RKF (Runge–Kutta–Fehlber​g), you need to transform your second order ODE into a first order one as:

$$\bar{u}^{'} = \bar{v}$$

$$\mu \bar{v}^{'} = \bar{p}^{'} + \rho \bar{v}$$

With boundary conditions: $$\bar{u}|_{y=-1} = 0$$ and $$\bar{u}|_{y=1} = 1$$.

In order to use shooting method you need to convert this boundary value problem (BVP) into an initial value problem (IVP) by replacing the second boundary condition ($$\bar{u}|_{y=1} = 1$$) with this one: $$\bar{u}^{'}|_{y=-1} = a$$, where $$a$$ is a unknown constant parameter here. Now, let's define this function:

$$F(a) = \bar{u}_{1}|_{y=1}(a) - \bar{u}|_{y=1}$$

Where $$\bar{u}_{1}|_{y=1}(a)$$ is the value of function obtained by solving IVP at $$y=1$$ and $$\bar{u}|_{y=1}$$ is the boundary condition of BVP above. If you find the root of $$F(a)$$, you will end up to having to finding that specific $$a$$ which would satisfy the second boundary condition, but that's what you want here as a solution for your BVP that started with it. In fact you solve the IVP by using RKF and then will find the root of $$F(a)$$ by some standard root finding methods such as Newton-Raphson or Bisection.

Update:

It's a Python implementation:

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import bisect
import matplotlib.pyplot as plt

mu = 1
rho = 1
u_b = 1

def p_prime(y):
return y

# function that returns dz/dt
def model(w,y):
dvdt = p_prime(y) / mu + (rho / mu) * w[1]
dudt = w[1]
dwdt = [dudt,dvdt]
return dwdt

def problem(a):

# initial condition
w0 = [0,a]

# time points
y = np.linspace(-1,1,100)

# solve ODE
w = odeint(model,w0,y)

return w

a = np.linspace(0,1,10)

y = np.linspace(-1,1,100)

for av in a:

w = problem(av)

# plot results
plt.plot(y,w[:,0],label=r'u'+' ,a = '+str(av))
#plt.plot(y,w[:,1],'r--',label=r'v')
plt.ylabel('u')
plt.xlabel('y')
plt.legend(loc='best')
plt.show()

def F(a):

w = problem(a)

return (w[y.shape[0]-1,0] - u_b)

a = np.linspace(0,1,100)

Fs = []

for av in a:
Fs.append(F(av))

# plot results
plt.plot(a,Fs)
plt.ylabel('F(a)')
plt.xlabel('a')
plt.show()

a_star = bisect(F,0,1)

print "The root of F(a): " + str(bisect(F,0,1))

w = problem(a_star)

# plot results
plt.plot(y,w[:,0],'b-',label=r'u'+' ,a = '+str(a_star))
plt.plot(y,w[:,1],'r--',label=r'v'+' ,a = '+str(a_star))
plt.ylabel('w')
plt.xlabel('y')
plt.legend(loc='best')
plt.show()


and these are outputs for $$\bar{u}$$ at different $$a$$ values, the variation of $$F(a)$$ at various $$a$$ values, and finally the $$a^{*}$$ which is the root of $$F(a)$$ and final result based on $$a^{*}$$:

The root of F(a): 0.469552906694