计算科学 - 高度非线性方程的有限差分 - 森林中的风 - 吾爱随笔录

高度非线性方程的有限差分 - 森林中的风

计算科学有限差分非线性方程微分方程

2021-12-24 02:53:30

基于 Navier-Stokes 方程和一些参数化，高度为 H 的森林内的水平稳态风 $u(z)$

a {(\frac{d u}{d z})}^{2} + b \frac{d u}{d z} \frac{d^{2} u}{d z^{2}} + c u + d \frac{d u}{d z} + e u^{2} + f = 0, for 0 < z < H .

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

系数随高度而变化，并且最初是给定的（我们可以根据需要对它们进行微分和积分）。 $a,b,c,d,e,f$ $z$

在地面： ,（也许没有必要在导数上指定 BC）。 $u(z = 0) = 0$ $\frac{du}{dz}(z=0) = 0$

在树冠顶部：。 $u(z = H) = U_H$

我正在尝试使用有限差分方案为求解这个方程，如果有人可以帮助我，那就太好了： $u(z)$

有限差分甚至是解决此类问题的好方法吗？
如果我使用经典表达式重写方程 etc... 我得到像这样的平方项，但我不知道如何从那里做什么。 $\frac{du}{dz} = \frac{u_{n+1}-u_{n-1}}{2h}, \frac{d^2u}{dz^2} =$ $u_{i+1}u_{i-1}$
我不知道如何正确使用牛顿法或皮卡德法，有没有更好的方法来重写方程？这样的变量？ $v = \frac{du}{dz}$

那时，我什至不确定我是否遗漏了一些明显的东西，或者这是一个非常困难的问题，任何帮助将不胜感激。

2个回答

我不能在这里谈论 FD 方法的使用，但我看不出它有任何行不通的原因。
有这些平方项和丑陋很好，它只是表明你的控制方程是非线性的。
一种方法是说你有一个残差 R，表示如下：在稳定状态下，R = 0。您可以添加伪时间导数，，并在伪时间离散化该术语以进入稳态答案。即，其中受 CFL 条件约束。这对应于前向欧拉时间步长。 $R (u) = a {(\frac{d u}{d z})}^{2} + b \frac{d u}{d z} \frac{d^{2} u}{d z^{2}} + c u + d \frac{d u}{d z} + e u^{2} + f$ $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$ $\frac{du}{dt}$ $u^{n + 1} = u^{n} + Δ t R (u)$ $u^{n+1} = u^n + \Delta t R(u)$ $\Delta t$

如果您想应用牛顿法，我们从 BDF1 和开始，其中不等于 0，而我们希望等于 0。如果我们对 u_n 进行 R(u_{n+1}) 的泰勒级数展开，我们：其中残差的导数是雅可比矩阵。然后我们可以将左项带入 RHS 并得到： $u_n$ $R(u_n)$ $R(u_{n+1})$ $R(u_{n+1})$ $u_n$

R (u_{n + 1}) = R (u_{n}) + \frac{\partial R}{\partial u_{n}} Δ u = 0

$R(u_{n+1}) = R(u_n) + \frac{\partial R}{\partial u_n}\Delta u = 0$

\frac{\partial R}{\partial u_{n}} Δ u = - R (u_{n})

$\frac{\partial R}{\partial u_n}\Delta u = -R(u_n)$ 现在，通过求解这个线性系统，我们实际上不会一步一步求解非线性控制方程，因为线性泰勒级数展开不足以快速求解非线性方程，但我们可以多次重复此迭代以求解稳态状态问题，您可能会注意到这与牛顿方法的相似之处更多。希望这有帮助！

除了同意有限差分方法是解决此问题的可接受方法外，这并没有具体回答您的问题。

相反，我将提出一种替代方法，假设您的主要目标是获得方程的解，而不是用数值方法进行试验。特别是对于像这样的具有挑战性的问题，我相信如果可能的话，使用该领域专家编写的软件而不是尝试编写自己的软件几乎总是更好。

如果你碰巧可以使用 MATLAB，我相信这个问题可以通过bvp4c函数解决。以下是我为解决您的方程式而编写的一些原型代码。但是我应该提前警告您，我为bvp4c能够获得解决方案的系数任意选择了一些值。我从我的数值实验中知道，对于许多值的组合，解决方案要么难以获得，要么不存在。

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

如果你没有 MATLAB，Python SciPy有一个类似于 bvp4c 的函数，你可以试试。

其它你可能感兴趣的问题

上一篇对广义相对论中计算能量条件的Python代码存疑下一篇是否存在关于计算科学实验如何声称有效性的“框架”？计算科学的科学方法？