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