Note: This tutorial was generated from an IPython notebook that can be downloaded here.

QuickstartΒΆ

This notebook was made with the following version of FEMpy:

[1]:
import FEMpy
FEMpy.__version__
[1]:
'1.0'

FEMpy can be used to solve Poisson equations on 1D and 2D domains. For this example, let us consider the problem:

\[-\frac{{\rm d}}{{\rm d} x} \left( e^x \frac{{\rm d}}{{\rm d} x} u(x)\right) = -e^x [\cos x -2\sin x - x \cos x - x \sin x]; x \in [0, 1]\]
\[u(0) = 0, u(1) = \cos(1)\]

To start, we will need to import our necessary aditional packages.

[2]:
import numpy as np

Let us define our coefficient and source functions,

[3]:
coefficient_function = lambda x: np.exp(x)
source_function = lambda x: -np.exp(x) * (np.cos(x) - 2*np.sin(x) - x*np.cos(x) - x*np.sin(x))

And our boundary condition,

[4]:
def dirichlet_function(x):
    if x == 0:
        return 0
    elif x == 1:
        return np.cos(1)

Under the Galerkin formulation, we can write our differiential equation as

\[\int_a^b c u' v' {\rm d}x = \int_a^b f v {\rm d}x\]

where \(v(x)\) is a test function. We then can choose our test and trial basis functions such that

\[u_h \in \{\phi \in C[a,b] \mid \phi(x) \text{ linear on each } [x_n, x_{n+1}]; (n = 1,2,\dotsc,N)\}\]

and

\[\int_a^b c u_h' v_h' {\rm d}x = \int_a^b f v_h {\rm d}x\]

for any \(v_h \in U_h\) Thus, we can define our test and trial basis function basis functions by:

[5]:
basis = FEMpy.IntervalBasis1D('linear')

We can set up our mesh using a step size of \(h=1/4\)

[6]:
mesh = FEMpy.Interval1D(left=0, right=1, h=1/4, basis_type='linear')

The boundary conditions are defined by:

[7]:
bcs = FEMpy.BoundaryConditions(mesh, boundary_types=('dirichlet', 'dirichlet'), dirichlet_fun=dirichlet_function)

Finally, we can input our mesh, basis functions, and boundary conditions into our Poisson equation then call our solver.

[8]:
poisson_eq = FEMpy.Poisson1D(mesh, fe_test_basis=basis, fe_trial_basis=basis, boundary_conditions=bcs)
poisson_eq.solve(coefficient_function, source_function)
[8]:
array([0.        , 0.24411715, 0.44112525, 0.55036422, 0.54030231])