Topic 15.2: Elliptic Partial-Differential Equations (Maple)

Introduction Notes Theory HOWTO Examples Engineering Error Questions Matlab Maple

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
);