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

Чтобы понять, как работает 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)

numerical-methods numerical-integration fipy
person homocomputeris    schedule 08.05.2019    source источник

Ответы (1)


arrow_upward
1
arrow_downward

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

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 (изменить из комментария)

person wd15    schedule 08.05.2019
comment
Также обратите внимание на изменение dx с 1 / nx на 1. / nx. Это не имеет значения в Py3k, но имеет значение в Python 2.7. - person jeguyer; 09.05.2019
comment
В Python 3 у меня нет PySparse, и когда я использую те же решатели из подмодуля решателей SciPy, они не дают ни NaN, ни нулей, хотя проблема не кажется сложной. - person homocomputeris; 09.05.2019
comment
Получить установку Python 2.7 с помощью решателей PySparse и Trilinos несложно: ctcms.nist.gov/fipy/INSTALLATION.html#recommended-method - person jeguyer; 09.05.2019