Finite difference for a highly nonlinear equation - The wind within the forest

Question

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\left(\frac{du}{dz}\right)^2 + b\frac{du}{dz} \frac{d^2u}{dz^2} + cu +d\frac{du}{dz} + eu^2 + f= 0,\quad\text{for}\quad 0<z<H. $$

The coefficients $a,b,c,d,e,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$, $\frac{du}{dz}(z=0) = 0$ (maybe it is not necessary to specify BC on the derivative).

At canopy top: $u(z = H) = U_H$.

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.

Answer 1

1. I can't speak to the use of the FD approach here, but I don't see any reason it wouldn't work.
2. Its fine to have those square terms and ugliness, it just indicates that your governing equations are nonlinear.
3. One way to approach this to say you have a residual, R, that is expressed like so: $$R(u) = a\left(\frac{du}{dz}\right)^2 + b\frac{du}{dz}\frac{d^2u}{dz^2} + cu + d\frac{du}{dz} + eu^2+f$$ And at steady state, R = 0. You could add a pseudo time derivative, $\frac{du}{dt}$, and discretize this term in pseudo time to march to the steady state answer. i.e., $$u^{n+1} = u^n + \Delta t R(u)$$ where $\Delta t$ is constrained by the CFL condition. This corresponds to forward euler time stepping.

Should you want to apply the newton method, we begin with BDF1 and with $u_n$, for which $R(u_n)$ does not equal zero, and a $R(u_{n+1})$ which we want to be equal to zero. If we take the taylor series expansion of $R(u_{n+1})$ about $u_n$ we get: $$R(u_{n+1}) = R(u_n) + \frac{\partial R}{\partial u_n}\Delta u = 0$$ Where the derivative of the residual is the jacobian matrix. We can then bring the left term to the RHS and get: $$\frac{\partial R}{\partial u_n}\Delta u = -R(u_n)$$ Now, by solving this linear system, we will not actually solve the nonlinear governing equations in one step, because a linear taylor series expansion is not sufficient to solve the nonlinear equation that quickly, but we can repeat this iteration many times to solve the steady state problem, and you may notice this bears more that a passing resemblance to Newton's Method. Hope this was helpful!

Answer 2

This doesn't specifically answer your question other than to agree that a finite difference method is an acceptable approach to this problem.

Instead, I'm going to propose an alternate approach assuming your main objective is to get a solution to your equation rather than experiment with numerical methods. Particularly for challenging problems like this one, I believe it is almost always better, if possible, to use software written by specialists in this area rather than try to write your own.

If you happen to have access to MATLAB, I believe this problem may be solvable with the bvp4c function. Below is some prototype code I wrote to solve your equation. But I should warn you in advance that I arbitrarily picked some values for the coefficients that allowed bvp4c to obtain a solution. I know from my numerical experiments that solutions are either difficult to obtain or don't exist for many combinations of values.

function answers_5_14_2019
dydx=@(x,y) f(x,y);
Uh=3;
bcFunc=@(ya,yb) bc(ya,yb,Uh);
yinit=@(x)[x; 1];
L=10;
solinit=bvpinit(linspace(0,L,20),yinit);
sol=bvp4c(dydx,bcFunc,solinit);
figure; plot(sol.x,sol.y(1,:)); grid;
end

function dydx = f(x,y)
a = .1; b = 1; c = 0; d = 1e-3*x; e = 0; f = 0;
  dydx = [y(2);
          -1/(b*y(2))*(a*y(2)^2 + c*y(1) +d*y(2) + e*y(1)^2 + f)];
end
% -----------------------------------------------------------------------
function res = bc(ya, yb, Uh)
% Boundary conditions
  res = [ya(1); yb(1)-Uh];  
end

If you don't have MATLAB, Python SciPy has a function similar to bvp4c that you can try.