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.


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

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:

# plot results

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

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

enter image description here

enter image description here

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.