maxima-bounces at math.utexas.edu wrote on 04/07/2009 06:14:02 PM:
> I have just commited a new graphic object (mesh) for package draw. Here
> is the associated documentation:
Thanks--your code worked great.
For anybody that is interested, I appended the code I used in class
today. It does a
finite difference solution to the heat equation. Presumably, this code
shows that
stability matters. To run this code, you'll need draw.lisp from CVS.
Again, thanks for this contribution.
Barton
/*
We'll solve the heat equation Ft = Fxx using an explicit finite
difference method and an implicit finite difference method. The
boundary conditions are
F(0,t) = 0,
F(1,t) = 0,
F(x,0) = 25 sin(pi x).
Let nx and nt be the number of spatial and temporal and spatial knots,
respectively. Let's solve the equation for t in [0, t_final]. We begin
by assigning values to nt, nx, and t_final.
*/
load("draw.lisp");
nt : 75$
nx : 75$
t_final : 0.3d0$
/*
The spatial and temporal knots are generated by the functions x_knot and
t_knot, respectively. These functions are
*/
x_knot( i ) := float(i / nx);
t_knot( i ) := float(t_final * i / nt);
/*
We'll save the solution in an array sol. We need to use the boundary
conditions to initialize the array sol.
*/
array(sol, nx+1,nt+1)$
for i : 0 thru nt do (
sol[0,i] : 0.0d0,
sol[nx,i] : 0.0d0)$
for i : 1 thru nx-1 do sol[i,0] : 25.0d0 * sin(float(%pi) * x_knot( i ))$
/*
Let dx and dt be the spatial and temporal knot spacings, respectively.
Thus
*/
dx : 1.0d0 / nx$
dt : t_final / nt$
/*
Let g = dt / dx^2. For stability, we need g < 1.
*/
g : dt / (dx * dx);
/*
Yikes! g = 22.5, so we expect the errors in the solution to grow
exponentially.
*/
a : 1 - 2 * g$
/*
We now fill in the array and plot the result
*/
for j : 1 thru nt do (
for i : 1 thru nx-1 do (
m : j - 1,
sol[i, j] : g * (sol[i - 1, m] + sol[i + 1, m]) + a*sol[i, m]))$
/*
For t in [0, final / 10], the solution looks OK:
*/
draw3d(
color = blue,
mesh(genmatrix(sol, nx, floor(nt/10),1,1),0,0,1, t_final/10),
xlabel = "x",
ylabel = "y",
surface_hide = true)$
/*
But for t in [0, t_final / 8], the graph starts to wiggle in a way that
doesn't match our expectation (on its own,
heat "spreads," it doesn't concentrate.)
*/
draw3d(
color = blue,
mesh(genmatrix(sol, nx, floor(nt/8),1,1),0,0,1, t_final/8),
xlabel = "x",
ylabel = "y",
surface_hide = true)$
/*
Let's try the implicit method. First, we'll set up the coefficient matrix
and compute its LU factorization.
*/
array(mm,nx,nx)$
for i : 0 thru nt do (
sol[0,i] : 0.0d0,
sol[nx,i] : 0.0d0)$
for i : 1 thru nx-1 do (
sol[i,0] : 25.0d0 * sin(float(%pi) * x_knot( i )))$
for i : 1 thru nx do ( for j : 1 thru nx do (
mm[i, j ] : 0.0d0))$
a : 1 + 2*g;
for i : 1 thru nx do mm[i , i] : a$
for i : 2 thru nx do (
mm[i,i-1] : -g,
mm[i-1,i] : -g)$
mmm : genmatrix(mm, nx,nx)$
mmm : lu_factor(mmm, 'floatfield)$
array(b, nx,1)$
for j : 1 thru nt do (
for i : 1 thru nx do (
b[ i,1 ] : sol[i,j-1]),
b[1,1] : b[1,1] + g * sol[0,j],
b[nx,1] : b[nx,1] + g * sol[nx,j],
bb : genmatrix(b, nx,1),
bb : lu_backsub(mmm, bb),
for i : 1 thru nx do (
sol[i, j] : bb[ i, 1]))$
/*
We win! For a fixed space point, the temperature steadily declines with
increasing time. This matches our
expectations. The moral is: more work better results (well, not always :)
*/
draw3d(
color = blue,
mesh(genmatrix(sol, nx, nt,1,1),0,0,1, t_final),
xlabel = "x",
ylabel = "y",
surface_hide = true)$