时间戳

首页 > 解决方案 > 如何在 FiPy 中求解 Euler-Bernoulli 梁方程？

numerical-methods - 如何在 FiPy 中求解 Euler-Bernoulli 梁方程？

问题描述

要了解 FiPy 的工作原理，我想求解具有固定端点的Euler–Bernoulli 梁方程：

w''''(x) = q(x,t),    w(0) = w(1) = 0,  w'(0) = w'(1) = 0.

为简单起见，让q(x,t) = sin(x).

如何在 FiPy 中定义和解决它？如何sin(x)针对方程中唯一的自变量指定源项？

from fipy import CellVariable, Grid1D, DiffusionTerm, ExplicitDiffusionTerm
from fipy.tools import numerix

nx = 50
dx = 1/nx
mesh = Grid1D(nx=nx, dx=dx)

w = CellVariable(name="deformation",mesh=mesh,value=0.0)

valueLeft = 0.0
valueRight = 0.0

w.constrain(valueLeft, mesh.facesLeft)
w.constrain(valueRight, mesh.facesRight)
w.faceGrad.constrain(valueLeft, mesh.facesLeft)
w.faceGrad.constrain(valueRight, mesh.facesRight)

# does not work:
eqX = DiffusionTerm((1.0, 1.0)) == numerix.sin(x)
eqX.solve(var=w)

标签： numerical-methodsnumerical-integrationfipy

解决方案

这似乎是您的问题的工作版本

from fipy import CellVariable, Grid1D, DiffusionTerm
from fipy.tools import numerix
from fipy.solvers.pysparse.linearPCGSolver import LinearPCGSolver
from fipy import Viewer
import numpy as np


L = 1.
nx = 500
dx = L / nx
mesh = Grid1D(nx=nx, dx=dx)

w = CellVariable(name="deformation",mesh=mesh,value=0.0)

valueLeft = 0.0
valueRight = 0.0

w.constrain(valueLeft, mesh.facesLeft)
w.constrain(valueRight, mesh.facesRight)
w.faceGrad.constrain(valueLeft, mesh.facesLeft)
w.faceGrad.constrain(valueRight, mesh.facesRight)

x = mesh.x
k_0 = 0
k_1 = -1
k_2 = 2 + np.cos(L) - 3 * np.sin(L)
k_3 = -1 + 2 * np.sin(L) - np.cos(L)
w_analytical = numerix.sin(x) + k_3 * x**3 + k_2 * x**2 + k_1 * x + k_0
w_analytical.name = 'analytical'

# does not work:
eqX = DiffusionTerm((1.0, 1.0)) == numerix.sin(x)
eqX.solve(var=w, solver=LinearPCGSolver(iterations=20000))
Viewer([w_analytical, w]).plot()
raw_input('stopped')

运行此之后，FiPy 解决方案似乎与分析结果非常接近。

OP 实施的两个重要变化。

  • 使用mesh.xwhich 是引用空间变量以在 FiPy 方程中使用的正确方法。

  • 指定求解器和迭代次数。该问题似乎收敛缓慢，因此需要大量迭代。根据我的经验，四阶空间方程通常需要良好的预处理器才能快速收敛。您可以尝试将 Trilinos 求解器包与 Fipy 一起使用，以使这项工作更好，因为它具有更广泛的可用预处理器。

  • 在 Python 2.7 中使用显式浮点数L来避免整数数学（从评论编辑）

推荐阅读