Hi,
I wrote a simple test script for a 1D ODE. I'm having trouble getting
the fipy solution to match the analytical solution. I'm pretty sure
the analytical solution is correct (but I guess I could be wrong). One
thing that really bothers me is that changing the number of points (nx
in the script below) drastically changes the solution fipy returns.
If you care about the physics of the problem: This a linearized beam
bending problem with the position and slope fixed at one end (x = 0)
and a dimensionless force 'f' applied at x = L. I was planning to add
a time dependent term in later, but I wanted to make sure I could
match a simpler problem with its analytical solution.
Thanks in advance for your help.
-Tony
Specs:
OS X 10.5.2
FiPy trunk
Numpy 1.0.4
Python 2.5.1
#===========
from __future__ import division
from fipy import *
Lx = 1.
nx = 100
dx = Lx / nx
mesh = Grid1D(dx=dx, nx=nx)
x = mesh.getCellCenters()[0]
k = 5
f = 5
y_analytical = CellVariable(mesh=mesh, name='analytical',
value = f / k**3 * (sin(k)*(cos(k*x) - 1) + cos(k)*(-sin(k*x) +
k*x)))
# fixed at x = 0; forced at x = Lx
BCs = (FixedValue(mesh.getFacesLeft(), value=0.),
FixedFlux(mesh.getFacesLeft(), value=0.),
NthOrderBoundaryCondition(mesh.getFacesRight(),
value=0., order=2),
NthOrderBoundaryCondition(mesh.getFacesRight(), value=f,
order=3))
# y_xxxx + k^2 * y_xx = 0
ElasticityTerm = ImplicitDiffusionTerm(coeff=(1.0, 1.0))
InertiaTerm = ImplicitDiffusionTerm(coeff=k**2)
eq = ElasticityTerm + InertiaTerm == 0.
solver = LinearLUSolver()
y = CellVariable(mesh=mesh, name='fipy', value=0.)
viewer = viewers.make(vars=(y, y_analytical))
eq.solve(var=y, boundaryConditions=BCs, solver=solver)
viewer.plot('./test')
