[MATLAB]: дискретное уравнение Пуассона

[DLL]

День добрый. Мне необходимо решить задачу Неймана для уравнение Пуассона при специфической дискретизации в квадрате.
Теоретически это можно написать и самому, но было бы проще воспользоваться готовым скриптом, например для стандартной дискретизации "крест".
Может кто-то видел учебный готовый скрипт в интернете? Немного пробежался гуглом, сходу не нашёл!

[vasya321]

DLL, вот скрипт решения ур. Пуассона ($\nabla^2 \varphi = -\rho/\varepsilon_0$) с нулевыми граничными условиями в квадратной области:

function phi = Poisson2D(H,dx,f_rho,eps)

% H - длина стороны квадрата
% dx - шаг сетки
% f_rho(x,y) - функция распределения плотности объемного заряда
% eps - заданная точность расчета

eps0 = 8.85E-12;

Nxe = H/dx + 2;
phi = zeros(Nxe);
phi_old = zeros(Nxe);
[X Y] = meshgrid(-dx/2:dx:H+dx/2);
rho = f_rho(X,Y);

iter = 1;
del = 1;
while del>eps
    for J=2:Nxe-1
        for I=2:Nxe-1
            phi(J,I) = 0.25*(phi_old(J+1,I)+phi_old(J,I+1)+phi_old(J-1,I)+phi_old(J,I-1)+rho(J,I)*dx^2/eps0);
        end
    end
    del = (phi(:)-phi_old(:))/phi(:);
    del = max(abs(del(~isnan(del))));
    fprintf('iter=%d, del=%G\n',iter,del);
    phi_old = phi;
    iter = iter+1;
end

Подойдет для начала?

[DLL]

Подойдет, спасибо!
Вопрос сразу: а зачем систему решать простой итерацией?
У меня MATLAB обращает матрицу для phi, быстрее чем простая итерация сходится

[vasya321]

Если Вы имеете в виду, почему не решить получившуюся после дискретизации исходного уравнения систему $A \varphi = B$ как $\varphi = A^{-1} B$, то никто не запрещает так делать и это действительно будет быстрее. Преимуществами решения итерационным методом являются:
1) используется меньше оперативной памяти, поскольку не хранятся нулевые элементы матрицы $A$,
2) не надо использовать разреженные матрицы,
3) более наглядное представление разностного шаблона и задание граничных условий, чем при заполнении матрицы $A$ и вектора $B$,
4) привычка решать данное уравнение именно этим методом.

Для ускорение работы функции выше можно использовать метод Гаусса-Зейделя, т.е. использовать формулу для phi(J,I):

phi(J,I) = 0.25*(phi_old(J+1,I)+phi_old(J,I+1)+phi(J-1,I)+phi(J,I-1)+rho(J,I)*dx^2/eps0);

или метод верхней релаксации, т.е. рассчитывать phi(J,I) как:

phi_tmp = 0.25*(phi_old(J+1,I)+phi_old(J,I+1)+phi(J-1,I)+phi(J,I-1)+rho(J,I)*dx^2/eps0);
phi(J,I) = omega*phi_tmp + (1-omega)*phi_old(J,I);

где $\mathrm{omega} \approx 1,85$.