Численное решение задачи Дирехле средствами matlab

Листинг pde/numgridborder.m

function B = regionborder(G)

% G - матрица, возвращаемая функцией numgrid

% B - то, что останется от G, если убрать внутренние точки

[m, n] = size(G);

inner = find(G);

border = []; % номера граничных точек

for k = [−1, 1, −m, m, −1−m, −1+m, 1−m, 1+m]

% точка - граничная,

% если по крайней мере одна из восьми соседних

% не принадлежит области

% for k = [-1, 1, -m, m]

% точка - граничная,

% если справа или слева или снизу или сверху есть точка,

% не принадлежащая области

border = [border; inner(G(inner + k) == 0)];

end;

B = zeros(size(G));

B(border) = G(border);

Решение задачи Дирихле реализует следующая функция.

Листинг pde/dirichlet.m

Задача Дирихле для уравнения Пуассона ∆u = fn(x, y), u|∂Γ = fnborder

function [X, Y, U ] = dirichlet(h, Region, fn, fnborder, varargin)

if nargin < 1 | | isempty(h)

h = 0.02;

end;

if nargin < 2 | | isempty(Region)

Region = ’L’;

end;

if nargin < 3 | | isempty(fn)

fn = @fun;

end;

if nargin < 4 | | isempty(fnborder)

fnborder = @funborder;

end;

x0 = 0;

x1 = 1;

n = (x1 - x0)/h + 3;

[X, Y ] = meshgrid(linspace(x0 − h, x1 + h, n));

U = zeros(n);

U(:) = NaN;

R = numgrid(Region, n); % нумерация узлов области

B = numgridborder(R); % узлы, лежащие на границе

Region = R > 0; % шаблон из нулей и единиц для все узлов области

Border = B > 0; % шаблон из нулей и единиц для узлов на границе

region = R(Region); % номера всех узлов в определенном порядке

border = R(Border); % номера узлов, лежащих на границе

nregion = length(region);

nborder = length(border);

A = -delsq(R)/h^2;

% правая часть уравнения

f = zeros(nregion, 1);

f (region) = feval(fn, X (Region), Y (Region), varargin{:});

% граничные условия

A(border,:) = sparse(1:nborder, border, 1, nborder, nregion);

f (border) = feval(fnborder, X (Border), Y (Border), varargin{:});

u = A\f;

U(Region) = u(region);

X = X(2: n - 1, 2: n - 1);

Y = Y (2: n - 1, 2: n - 1);

U = U(2: n - 1, 2: n - 1);

Clf

Shg

surfc(X, Y, U);

shading interp;

View(-63, 17)

function f = fun(x, y)

f = 50*(4x-2y)

function f = funborder(x, y)

f = 0;

% f = 2x-3y;

% f = sin(pi*x) - sin(pi*y);

Приложение 7