Fipy/Simple Transient Problem

This page covers the solution of a simple ODE using Fipy.

The Equation

This script solves the following differential equation:

$ \frac{dy}{dt} = -1.8 y - 1.0 $

with initial condition

$ y(t=0) = 1.0 $

which has the exact solution:

$ y(t) = 1.5556 \exp( -1.8 t ) - 0.5556 $

The Script

Import Statements

We begin with our usual import statements, which allow us to use Fipy for solving the equation and Matplotlib for plotting results.

from numpy import *
import matplotlib.pylab as plt
from fipy import *

# Set up figure for plotting
fig = plt.figure()
ax = fig.add_subplot(111)

Header

Here's a header, to document what this file does.

"""
Simple Transient Problem

Charles Reid
January 2014

This script solves the ODE:

dy/dt = - 1.8*y - 1.0

y(t=0) = 1.0

This has analytical solution:

y(t) = 1.5556 exp(-1.8 t) - 0.5556

(see http://wolfr.am/1eveai5)
"""

Input Parameters

The first thing we'll do is set some parameters, which we will then use to construct Fipy objects:

##############################
# Set Fipy parameters

# number of cells in mesh
nx = 1

# timesteps
deltat = 0.1
nt = 100

# initial value
iv = 1.0

Fipy Objects

Now assemble Fipy objects using user parameters

##############################
# Assemble Fipy objects

# create the mesh
mesh = Grid1D(nx=nx,dx=1.0)

# create the variable
y = CellVariable(mesh=mesh, value=iv, hasOld=1, name='y')

# create the equation
RHS = -1.8*y - 1.0
eq = TransientTerm() == RHS

# solver
from fipy.solvers import LinearPCGSolver
solver = LinearPCGSolver()
#from fipy.solvers import DefaultAsymmetricSolver
#solver = DefaultAsymmetricSolver()

# timesteps and steps 
dt = Variable(deltat)
steps = nt
t = dt.value * steps
solutions = zeros([steps+1,nx])
solutions[0,:] = y.value

Solution Loop

Loop over each step and solve the differential equation over that step:

# loop over timesteps 
for step in [i+1 for i in range(steps)]:
    y.updateOld()
    eq.solve(y,
            solver = solver,
            dt=dt)

    # store solutions
    solutions[step] = y.value

Plot Computed and Exact Solution

ts = linspace(0,t,steps)

# plot computed solution
ax.plot(ts,solutions[:step],'bo',label='Computed')#arange(len(y.value)),y.value)

# plot exact solution
y2 = 1.5556*exp(-1.8*ts) - 0.5556 
plt.hold(True)
ax.plot(ts,y2,'r-',label='Exact')

ax.legend(); plt.show(); plt.draw()

The Result

This should result in the following plot: