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

The Equation

This script solves the following differential equation:

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

with initial condition

$y(t=0)=1.0$

which has the exact solution:

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

More information about this problem is available via Wolfram Alpha.

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)

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:

The Py File

http://files.charlesmartinreid.com/FipySimpleTransient.py

Fipy/Simple Transient Problem

From charlesreid1

Contents