"""Solve u'=-a*u, u(0)=I, for t in (0,T] with steps of dt."""

dt = float(dt) # avoid integer division

Nt = int(round(T/dt)) # no of time intervals

T = Nt*dt # adjust T to fit time step dt

u = zeros(Nt+1) # array of u[n] values

t = linspace(0, T, Nt+1) # time mesh

u[0] = I # assign initial condition

for n in range(0, Nt): # n=0,1,...,Nt-1

u[n+1] = (1 - (1-theta)*a*dt)/(1 + theta*dt*a)*u[n]

"""

Run a case with the solver, compute error measure,

and plot the numerical and exact solutions (if makeplot=True).

"""

u, t = solver(I, a, T, dt, theta) # Numerical solution

u_e = exact_solution(t, I, a)

e = u_e - u

E = sqrt(dt*sum(e**2))

if makeplot:

figure() # create new plot

t_e = linspace(0, T, 1001) # very fine mesh for u_e

u_e = exact_solution(t_e, I, a)

plot(t, u, 'r--o') # red dashes w/circles

plot(t_e, u_e, 'b-') # blue line for u_e

legend(['numerical', 'exact'])

xlabel('t')

ylabel('u')

title('Method: theta-rule, theta=%g, dt=%g' % (theta, dt))

theta2name = {0: 'FE', 1: 'BE', 0.5: 'CN'}

savefig('%s_%g.png' % (theta2name[theta], dt))

show()

import argparse

parser = argparse.ArgumentParser()

parser.add_argument('--I', '--initial_condition', type=float,

default=1.0, help='initial condition, u(0)',

metavar='I')

parser.add_argument('--a', type=float,

default=1.0, help='coefficient in ODE',

metavar='a')

parser.add_argument('--T', '--stop_time', type=float,

default=1.0, help='end time of simulation',

metavar='T')

parser.add_argument('--makeplot', action='store_true',

help='display plot or not')

parser.add_argument('--dt', '--time_step_values', type=float,

default=[1.0], help='time step values',

metavar='dt', nargs='+', dest='dt_values')

parser = define_command_line_options()

args = parser.parse_args()

print 'I={}, a={}, T={}, makeplot={}, dt_values={}'.format(

args.I, args.a, args.T, args.makeplot, args.dt_values)

I, a, T, makeplot, dt_values = read_command_line()

for theta in 0, 0.5, 1:

for dt in dt_values:

E = explore(I, a, T, dt, theta, makeplot)