pde_laplace.mw

FDM fuer die Laplace-Gleichung (Dirichlet-Problem) 

> restart:
 

> # Solving Laplace's equation on [0, 1] × [0, 1] where the boundary value is
# defined by a function f(x, y) may be done as follows:
 

> N :=32:                          # the number of intervals
 

> f := (x, y) -> sin(4.0*x*y):
 

> # assign the boundary values on [0, 1] x [0, 1]
 

> for i from 0 to N do
   du[i,0] := f(i/N, 0);
   du[0,i] := f(0, i/N);
   du[i,N] := f(i/N, 1);
   du[N,i] := f(1, i/N);
end do:
 

> M := Matrix( (N - 1)*(N - 1), (N - 1)*(N - 1), datatype=float[8] ):
b := Vector( (N - 1)*(N - 1), datatype=float[8] ):
 

> for i from 1 to N - 1 do
   for j from 1 to N - 1 do
       posn := (N - 1)*(i - 1) + j;
       M[posn, posn] := -4;

       for k from -1 to 1 by 2 do            if assigned( du[i + k, j] ) then
               b[posn] := b[posn] - du[i + k, j];
           else
               M[posn, posn + (N - 1)*k] := 1;
           end if;
       end do;

       for k from -1 to 1 by 2 do
           if assigned( du[i, j + k] ) then
              b[posn] := b[posn] - du[i, j + k];
           else
              M[posn, posn + k] := 1;
          end if;
       end do;
   end do;
end do:

u := LinearAlgebra:-LinearSolve( M, b ):   # solve the system of linear equations
 

> # plot the result
plots[display](
   plots[pointplot3d](
       [seq( seq( [i/N, j/N, u[(i - 1)*(N - 1) + j]], i = 1..N - 1 ), j = 1..N - 1 )],
       symbol = circle, color = black
   ),
   # the boundary values    
   plots[spacecurve]( [0, t, f(0, t)], t = 0..1, thickness = 2, colour = red ),
   plots[spacecurve]( [t, 0, f(t, 0)], t = 0..1, thickness = 2, colour = red ),
   plots[spacecurve]( [1, t, f(1, t)], t = 0..1, thickness = 2, colour = red ),
   plots[spacecurve]( [t, 1, f(t, 1)], t = 0..1, thickness = 2, colour = red ),
   axes = framed, axis[1,2]=[mode=linear]
);
 

Plot 

>