Как решить уравнение пучка Эйлера – Бернулли в FiPy? - PullRequest
Новая
       74

Как решить уравнение пучка Эйлера – Бернулли в FiPy?

1 голос
/ 08 мая 2019

Чтобы понять, как работает FiPy, я хочу решить уравнение пучка Эйлера – Бернулли с фиксированными конечными точками: 

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)

1 Ответ

1 голос
/ 08 мая 2019

Вот то, что кажется рабочей версией вашей проблемы 

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.x, который является правильным способом ссылки на пространственную переменную для использования в уравнениях FiPy.

  • Указание решателя и количества итераций.Проблема кажется медленной, чтобы сходиться, поэтому потребовалось много итераций.По моему опыту, пространственные уравнения четвертого порядка часто нуждаются в хороших предварительных кондиционерах, чтобы быстро сходиться.Вы можете попробовать использовать пакет решателя Trilinos с Fipy, чтобы улучшить эту работу, так как он имеет более широкий диапазон доступных предварительных кондиционеров.

  • Используйте явное значение с плавающей запятой для L, чтобы избежать целочисленной математики вPython 2.7 (редактировать из комментария)

...