

%
% Solve the Dirichlet problem for Laplace's equation using
% Leibman's five-point finit-differenc approximation. 
% The goal ?- go1 is a normal example, while the goal ?- go2
% shows output constraints for a small region where the boundary conditions
% are not specified.
%
:- use_package(clpq).
:- use_module(library(format)).

laplace([_, _]).
laplace([H1, H2, H3|T]):-
	laplace_vec(H1, H2, H3), 
	laplace([H2, H3|T]).

laplace_vec([_, _], [_, _], [_, _]).
laplace_vec([_TL, T, TR|T1], [ML, M, MR|T2], [_BL, B, BR|T3]):-
	B + T + ML + MR - 4 * M .=. 0, 
	laplace_vec([T, TR|T1], [M, MR|T2], [B, BR|T3]).

printmat([]).
printmat([H|T]):-
	printvec(H), 
	printmat(T).

printvec([]):- nl.
printvec([H|T]):-
	printrat(H), 
	printvec(T).

printrat(rat(N,D)) :- !,
        X is N/D,
        format(" ~2f",X).
printrat(N) :-
        X is N*100,
        format(" ~2d",X).

go1:-
	X =  [
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
    [100, _, _, _, _, _, _, _, _, _, 100], 
    [100, _, _, _, _, _, _, _, _, _, 100], 
    [100, _, _, _, _, _, _, _, _, _, 100], 
    [100, _, _, _, _, _, _, _, _, _, 100], 
    [100, _, _, _, _, _, _, _, _, _, 100], 
    [100, _, _, _, _, _, _, _, _, _, 100], 
    [100, _, _, _, _, _, _, _, _, _, 100], 
    [100, _, _, _, _, _, _, _, _, _, 100], 
    [100, _, _, _, _, _, _, _, _, _, 100], 
    [100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100]
    ], 
	laplace(X),
	printmat(X).

% Answer:
% 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
% 100.00 51.11 32.52 24.56 21.11 20.12 21.11 24.56 32.52 51.11 100.00
% 100.00 71.91 54.41 44.63 39.74 38.26 39.74 44.63 54.41 71.91 100.00
% 100.00 82.12 68.59 59.80 54.97 53.44 54.97 59.80 68.59 82.12 100.00
% 100.00 87.97 78.03 71.00 66.90 65.56 66.90 71.00 78.03 87.97 100.00
% 100.00 91.71 84.58 79.28 76.07 75.00 76.07 79.28 84.58 91.71 100.00
% 100.00 94.30 89.29 85.47 83.10 82.30 83.10 85.47 89.29 94.30 100.00
% 100.00 96.20 92.82 90.20 88.56 88.00 88.56 90.20 92.82 96.20 100.00
% 100.00 97.67 95.59 93.96 92.93 92.58 92.93 93.96 95.59 97.67 100.00
% 100.00 98.89 97.90 97.12 96.63 96.46 96.63 97.12 97.90 98.89 100.00
% 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00

go2([B31, M32, M33, B34, B42, B43, B12, B13, B21, M22, M23, B24]) :-
	laplace([
    	[_B11, B12, B13, _B14], 
    	[B21, M22, M23, B24], 
    	[B31, M32, M33, B34], 
    	[_B41, B42, B43, _B44]
    ]).
    	  

% Answer:
% 
% B34.=. -4*M22+B12+B21+4*M33-B43,
% M23.=.4*M22-M32-B12-B21,
% B31.=. -M22+4*M32-M33-B42,
% B24.=.15*M22-4*M32-4*B12-4*B21-M33-B13 ? 
