# Shooting method - Matlab ODE

I'm trying to solve these equations of hypersonic adiabatic flow over a flat plate. I did all the simplifications and got these equations for the stagnation point flow. $$\left(Cf''\right)' + f f'' = \left(f'\right)^2-g$$ $$\left(\frac{C}{p_r} g'\right) + fg'=0$$

For a calorically perfect gas, $g=h/he$, and $he = \rm stagnation ~enthalpy$

$p_r$ is a constant and $C=g^{-1/3}$

$f(0)=0$ (over the streamline of the body)

$f'(0)=0$ (this is $u$, the velocity over the body which is null)

$f''(0)=?$ white suggests to try as a first value $0.664/\sqrt{C_w}$

I'm doing this for the adiabatic case so

$g(0)=g_{wall}$

$g'(0)=0$ (no thermic flux)

UPDATE

Control conditions:
$f'(\delta)=1$ and $g(\delta)=1$ I want to integrate from $\delta=0$ to $\delta=6$, in the end, if $f'(\delta)$ is different than $1$ I'll have to change the value assumed for $f''(0)$.

How can I solve these equations, someone suggested me using the shooting method but how could I code this in matlab? I would like just a help on how to implement the equations on matlab, not a final solution.

• Similar kind of question is found in scicomp.stackexchange.com/questions/21498/…. Here you should assume different values of $f''(0)$ to get $g'(0)=0$. Still having problem pls. comment here.
– AGN
Jan 4 '16 at 1:15
• Are there any wrong/missing boundary conditions? You can't apply the shooting method unless you have a boundary condition at a location not at x=0. The whole premise of the shooting method is you treat a boundary value problem like an initial value problem and "shoot" from the one boundary, say the left, using the left conditions as initial conditions, and then compare the results to whatever the right boundary conditions are. If the results don't match up with the right BCs, you modify the unknown left BCs until you get there. Jan 4 '16 at 2:19
• @ArunGovindNeelanA Yes that's exactly what I need todo, different values of f''(0) to get g'(0)=0. But I have 2 equations, since h is unknown. In the other question and materials you suggested there was only one function, I have 2, f and g Jan 4 '16 at 10:36
• Right now I'm slightly busy, I will provide some details later. I think you are clear in procedure and only problem is computation. First split the third order ode into three first order ode then you can use RK methods. Then feed four ode with IC to ode45, for that you shall refer www3.nd.edu/~nancy/Math20750/Demos/3dplots/dim3system.html
– AGN
Jan 5 '16 at 6:36
• I'm guessing $\delta$ is fixed? I'm not sure about specifying both $f'$ and $g$ at $\delta$; looks like too many boundary conditions. Also if you set $f''(0)=0$ then you wont have a shooting parameter to use. Jan 7 '16 at 16:49

To change this into a system of first order ODEs, you can define

$$Y(x;s) := \left( \begin{split} f''\\ f'\\ f\\ g'\\ g \end{split} \right) = \left(\begin{split} Y_1\\ Y_2\\ Y_3\\ Y_4\\ Y_5 \end{split} \right)\,.$$ Your first equation then gives $$CY_1' + C'Y_1 + Y_3Y_5= (Y_2)^2 - Y_5\,,$$ or, with $C=g^{-1/3}$, $C' = -\frac13g^{-4/3}g'$, $$Y_5^{-1/3}Y_1' - \tfrac13Y_5^{-4/3} Y_4Y_1 + Y_3Y_5= (Y_2)^2 - Y_5\,,$$ and the second differentiated is $$\left(\frac{Y_5^{-1/3}}{p_r}+Y_3\right)Y_4' + \left(\frac{-\frac13Y_4Y_5^{-4/3}}{p_r}+Y_2\right)Y_4 = 0\,,$$ with the additional equations $$Y_2' = Y_1, \quad Y_3' = Y_2, \quad Y_5' = Y_4\,.$$

Bringing these together, gives $$Y'(x;s) = \left(\begin{split} Y_5^{1/3}\left(\tfrac13Y_5^{-4/3} Y_4Y_1 - Y_3Y_5 + (Y_2)^2 - Y_5\right)\\ Y_1\\ Y_2\\ - \left(\frac{Y_5^{-1/3}}{p_r}+Y_3\right)^{-1} \left(\frac{-\frac13Y_4Y_5^{-4/3}}{p_r}+Y_2\right)Y_4\\ Y_4 \end{split}\right)$$ and initial conditions $$Y(0;s) := \left( \begin{split} s\\ 0\\ 0\\ 0\\ g_{\text{wall}} \end{split} \right) \,.$$

You can use an ODE solver to integrate this the find $Y(\delta)$ for some $\delta>0$, and change $s$ to match your additional condition $$F(s) := Y_2(\delta;s) - 1 = 0\,,$$ or $$F(s) := Y_5(\delta;s) - 1 = 0\,.$$

This could be by bisection (find 2 values $s_1$ and $s_2$ with $F(s_1)<0$, $F(s_2)>0$, and then home in by testing $\tilde{s}=0.5(s_1+s_2)$ and replace $s_1$ or $s_2$ as appropriate) or even Newton's method (more work is required to get $F_s(s)$, but it's possible).

• It might be cleaner to define $Y_4 = C$, not sure. Jan 6 '16 at 13:43
• My bad, I didn't give the control conditions. n=f'(delta)=1 and g(delta)=1 I want to integrate from n=0 to n=6, in the end, if f'(delta) is different than 1 I'll have to change the value assumed for f''(0) Jan 7 '16 at 15:12
• I'm still struggling with the second equation $(C/p_r +f)g'=0$. Is this correct? Could you do something like saying in general $g' \neq 0$ so $C/p_r+f = 0$? In which case you could differentiate it to get an expression for $g'$ in terms of $Y_2$ and eliminate it from the equation for $Y_1'$ and carry on with the shooting method from there. Jan 7 '16 at 16:22
• Maybe you could try to minimise both your conditions with $F(s) := (Y_2(\delta;s)-1)^2+(Y_5(\delta;s)-1)^2$, but this is a different problem Jan 8 '16 at 11:43
• Does this make sense now? Hopefully it's now in a form where you can use the previous answers. Note you have to have $\delta >0$. Jan 12 '16 at 10:11