|
|
|
from numpy import *
def solver(I, a, T, dt, theta):
"""Solve u'=-a*u, u(0)=I, for t in (0,T] with steps of dt."""
Nt = int(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]
return u, t
Note about the for loop: range(0, Nt, s) generates all integers
from 0 to Nt in steps of s (default 1), but not including Nt (!).
Sample call:
u, t = solver(I=1, a=2, T=8, dt=0.8, theta=1)
|
|
|