4534.sagews

This worksheet is based on the description in Wikipedia on the Crank-Nicolson Method.

# Crank-Nicholson used to solve the 1 dimensional linear diffusion equation:  u_t = p*u_xx   where p is constant.
#define a test function for the bottom. I think of this as an initial concentration (of salt in water, heat in a metal bar) function
#over time this concentration changes according to the u_t = p*u_xx.  The concentration at the ends is forced by the side conditions
#u(t,0)=lft(t) and u(t,1)=rt(t).   
def f(x):
    return(x-x*cos(float(pi*.5)*x))
    #return(sin(float(pi)*x)+x)
    #return x
def lft(t):
    return(t)
def rt(t):
    return(1.0-t)    
def umat(f=x,lft=0,rt=1,N=16,M=4,p=1):
    dx=float(1/M);
    dt=float(1/N);
    r=p*.5*dt/dx^2;print(r)
    #set up the table of approximate values for u with the given values on the bottom and sides
    u=matrix(RR,N+1,M+1);
    for i in range(0,M+1):
        u[0,i]=f(i*dx)
    for i in range(1,N+1):
        u[i,0]=lft(i*dt)
    for i in range(1,N+1):
        u[i,M]=rt(i*dt)
    #initialize a,b,c and righthand side d
    a=vector(RR,{M:0})
    b=vector(RR,{M:0})
    c=vector(RR,{M:0})
    d=vector(RR,{M:0})
    #now compute the rows of u, using Crank-Nicholson equations for 1D diffusion and the tridiagonal system solution algorithm.     
    for k in range(1, N+1):
        #set up the coefficient matrix 
        for i in range(1,M):
            a[i]=-r
            b[i]=1+2*r
            c[i]=-r
        a[1]=0
        c[M-1]=0
        #set up the right hand side d
        for i in range(1, M): 
            d[i]=(r*u[k-1,i-1]+(1-2*r)*u[k-1,i]+r*u[k-1,i+1])
        d[M-1]=d[M-1]+r*u[k,M]
        d[1]=d[1]+r*u[k,0]    
        #the forward elimination step
        c[1]=c[1]/b[1]; d[1]=d[1]/b[1]
        for i in range(2,M-1):
            c[i]=c[i]/(b[i]-c[i-1]*a[i])
            d[i]=(d[i]-d[i-1]*a[i])/(b[i]-c[i-1]*a[i])
        d[M-1]=(d[M-1]-d[M-2]*a[M-1])/(b[M-1]-c[M-2]*a[M-1])
        #print(c);print(d);
        #the back substitution step
        u[k,M-1]=d[M-1]
        for i in range(0,M-2):
            u[k,M-2-i]=d[M-2-i]-c[M-2-i]*u[k,M-2-i+1]
    return(u)
    
# we can use plot_list to draw a picture of the concentration function at any time.  
def plotrow(u,f,i):
    if i > u.nrows():
        return('no row here')
    M=u.ncols()
    V=[]
    for j in range(0,M):
        V.append([j/(M-1),u[i,j]])
    return(list_plot(V)+plot(f,(x,0,1)))

# compute  an 11 by 101 table of approximate value of u.
w=umat(f,lft,rt,M=10,N=100,p=.5);#w

0.250000000000000

a = animate([plotrow(w,f,i) for i in range(0,100,5)], xmin=0, xmax=1, figsize=[2,1])

show(a)

#we can use list_plot3d and get an idea of what the distribution surface looks like
list_plot3d(w)

#We can also use point3d to draw piecewise linear curves thru the points in each row of the table

def plottable(u):
    M=u.ncols()
    N=u.nrows()
    L=list_plot3d(u)
    for r in range(0,N):
        V=[]
        for j in range(0,M):
            V.append((r,j,u[r,j]))
        L=L+line(V,color="red")
    return(L)

show(plottable(w))

animate?

File: /usr/local/lib/sage/sage-4.7.2/local/lib/python2.6/site-packages/sage/plot/animate.py

Type: <type ‘type’>

Definition: animate( [noargspec] )

Docstring:

Return an animation of a sequence of plots of objects.

INPUT:

• v - list of Sage objects. These should preferably be graphics objects, but if they aren’t then plot is called on them.
• xmin, xmax, ymin, ymax - the ranges of the x and y axes.
• **kwds - all additional inputs are passed onto the rendering command. E.g., use figsize to adjust the resolution and aspect ratio.

EXAMPLES:

sage: a = animate([sin(x + float(k)) for k in srange(0,2*pi,0.3)],
...                xmin=0, xmax=2*pi, figsize=[2,1])
sage: a
Animation with 21 frames
sage: a[:5]
Animation with 5 frames
sage: a.show()          # optional -- ImageMagick
sage: a[:5].show()      # optional -- ImageMagick

The show function takes arguments to specify the delay between frames (measured in hundredths of a second, default value 20) and the number of iterations (default value 0, which means to iterate forever). To iterate 4 times with half a second between each frame:

sage: a.show(delay=50, iterations=4) # optional -- ImageMagick

An animation of drawing a parabola:

sage: step = 0.1
sage: L = Graphics()
sage: v = []
sage: for i in srange(0,1,step):
...       L += line([(i,i^2),(i+step,(i+step)^2)], rgbcolor=(1,0,0), thickness=2)
...       v.append(L)
sage: a = animate(v, xmin=0, ymin=0)
sage: a.show() # optional -- ImageMagick
sage: show(L)

TESTS: This illustrates that ticket #2066 is fixed (setting axes ranges when an endpoint is 0):

sage: animate([plot(sin, -1,1)], xmin=0, ymin=0)._Animation__kwds['xmin']
0

We check that Trac #7981 is fixed:

sage: a = animate([plot(sin(x + float(k)), (0, 2*pi), ymin=-5, ymax=5)
...            for k in srange(0,2*pi,0.3)])
sage: a.show() # optional -- ImageMagick