5
$\begingroup$

Based on the Navier-Stokes equations and a few parameterizations, the horizontal steady-state wind $u(z)$ within a forest of height $H$ satisfies: $$ a\Big(\frac{du}{dz}\Big)^{\!2} + b\frac{du}{dz} \frac{d^2u}{dz^2} + cu +d\frac{du}{dz} + eu^2 + f= 0\:\text{ for }0<z<H. $$

The coefficients $a$ to $f$ vary with the altitude $z$ and are given initially (we can differentiable and integrate them as many times as needed).
At ground level: $$ u|_{z = 0} = 0, \quad\frac{du}{dz}\Big|_{z=0} = 0. $$ At canopy top: $$ u|_{z = H} = U_H, \quad\frac{du}{dz}\Big|_{z=H} = K\text{ (constant)} $$ I am trying to solve this equation for $u(z)$ using a finite difference scheme, it would be great if someone could help me:

  1. Are Finite Differences even a good approach for this kind of problem ?

  2. If I rewrite the equation using the classical expressions $\frac{du}{dz} = \frac{u_{n+1}-u_{n-1}}{2h}, \frac{d^2u}{dz^2} =$ etc... I obtain square terms like $u_{i+1}u_{i-1}$ and I do not how what to do from there.

  3. I do not know how to use the Newton method or the Picard method correctly, is there a better way to rewrite the equation ? Using variables like $v = \frac{du}{dz}$ for example ?

At that point, I am not even sure if I am missing something obvious or if this is a really hard problem, any help would be greatly appreciated.

Thanks

Improve this question
$\endgroup$
3
  • $\begingroup$ Are you sure that you want the boundary conditions for both the value and the derivative for a second order ODE? At $0$ you have a singularity, of course, but at $H$ the IVP is locally well-posed (if $K\ne 0$) and if the height is small enough and the coefficients are tame enough, you may arrive to $0$ with data incompatible with the boundary values there, so no numerical scheme will give you anything meaningful. $\endgroup$
    – fedja
    Commented May 15, 2019 at 2:22
  • $\begingroup$ You are probably right, there is no need to have BC on the derivatives. What do you mean by "if the height is small enough and the coefficients are tame enough, you may arrive to 0 with data incompatible " ? What is arriving to 0 ? $\endgroup$
    – Matt
    Commented May 15, 2019 at 6:27
  • $\begingroup$ "What is arriving to 0?" The solution of the initial value problem at $H$. Also, do you expect some a priori properties of the solution? (I would expect it at least to be increasing, perhaps even convex, but it may be too naive). Actually, if you could give a typical example of the coefficients you are dealing with, that might make things clearer too. $\endgroup$
    – fedja
    Commented May 15, 2019 at 9:37

1 Answer 1

Reset to default
1
$\begingroup$

As noted already in the comments, your boundary conditions seem off. Note that generically for a second-order BVP one expects to impose only two boundary conditions; you have 4.

Once you’re sure you’ve formulated the problem correctly, finite differences are a good first approach. A finite difference discretization will leave you with a system of nonlinear algebraic equations to solve, and I would suggest using a library for that (e.g. fsolve in Matlab or scipy). You could also directly use a library that solves BVPs.

To get a basic understanding of this kind of thing, I recommend reading the first two chapters of LeVeque's book on finite differences which includes an example of what can happen if you impose boundary conditions that make the problem ill-posed.

Since your problem is probably convection-dominated, you may find it necessary to use upwinded (one-sided) finite differences for the convective terms in order to avoid oscillations.

Finally, I recommend posting questions like this on scicomp.SE, as there are many experts that read that site but aren't on Mathoverflow.

Improve this answer
$\endgroup$
5
  • 1
    $\begingroup$ "I recommend posting questions like this on scicomp.SE" I would agree if not for the answers Matt got there now. Nobody there even bothered to ask him for the specifics of the problem or give a momentary thought to the existence and uniqueness issues yet, while it is rather clear that the problem is ill-posed at least in some cases. The advice he's getting there so far is the worst possible one: "just type the problem into some existing software or apply some generic scheme and hope for the best". I hope that will change, but right now I feel that even my stupid questions make more sense :( $\endgroup$
    – fedja
    Commented May 15, 2019 at 23:40
  • 1
    $\begingroup$ Respectfully, f you come asking how to compute solutions and you haven’t thought about well-posedness, I think you need to blame yourself and not the internet. And if your main concern is well-posedness then scicomp is probably not the best place. For computing solutions, existing software is often a good place to start; that doesn’t mean you should turn off your brain when you use it. $\endgroup$
    – David Ketcheson
    Commented May 16, 2019 at 4:08
  • $\begingroup$ Nevertheless, kudos to you @fedja for answering the question that should have been asked rather than the one that actually was asked. $\endgroup$
    – David Ketcheson
    Commented May 16, 2019 at 4:25
  • $\begingroup$ "if you come asking how to compute solutions and you haven’t thought about well-posedness, I think you need to blame yourself and not the internet". Of course, but when talking mathematics, the point is not to find the guilty person, but to figure out what the problem really is and how to solve it best, so "answering the question that should been asked" is sort of an industry standard in such cases. I do not see why scicomp should be any different. You shouldn't turn your brain off when following an advice, no doubt, but you should also try to turn it on as much as you can, when giving it :) $\endgroup$
    – fedja
    Commented May 16, 2019 at 5:42
  • $\begingroup$ Anyway, I disagree with you less than you may think. Maybe I just remember the times when I myself couldn't pose my questions in a good way a bit better, that's all. $\endgroup$
    – fedja
    Commented May 16, 2019 at 5:49

You must log in to answer this question.