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)=0$ (no thermic flux)


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.

  • 1
    $\begingroup$ 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. $\endgroup$ Jan 4, 2016 at 1:15
  • $\begingroup$ 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. $\endgroup$
    – spektr
    Jan 4, 2016 at 2:19
  • $\begingroup$ @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 $\endgroup$ Jan 4, 2016 at 10:36
  • $\begingroup$ 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 $\endgroup$ Jan 5, 2016 at 6:36
  • $\begingroup$ 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. $\endgroup$
    – Steve
    Jan 7, 2016 at 16:49

1 Answer 1


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

  • $\begingroup$ It might be cleaner to define $Y_4 = C$, not sure. $\endgroup$
    – Steve
    Jan 6, 2016 at 13:43
  • $\begingroup$ 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) $\endgroup$ Jan 7, 2016 at 15:12
  • 1
    $\begingroup$ 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. $\endgroup$
    – Steve
    Jan 7, 2016 at 16:22
  • $\begingroup$ 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 $\endgroup$
    – Steve
    Jan 8, 2016 at 11:43
  • $\begingroup$ 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$. $\endgroup$
    – Steve
    Jan 12, 2016 at 10:11

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.