Эллиптическое ДУ в частных производных

@tynec02 · Регистрация: 20.11.2013

Author24 — интернет-сервис помощи студентам

Здравствуйте, уважаемые форумчане. Буду благодарен за помощь со следующей задачей.

Собственно, условие:

В сети нарыл замечательную книжку, в которой нашел код для Решения двухмерного уравнения Пуассона методом конечных разностей, что, я думаю, мне и нужно. После адаптации к задаче, он выглядел так:

Matlab M

% u_xx + u_yy = 0, x,y [0,1], V1=xy, V2=x+y
 
clc
clear
close
 
x0=0; % x0 - начальная координата области решения по оси х;
xn=1; % xn - конечная координата области решения по оси х;
n=50; % n - число точек координатной сетки вдоль оси х; 
y0=0; % y0 - начальная координата области решения по оси y; 
ym=1; % ym - конечная координата области решения по оси y; 
m=30; % m - число точек координатной сетки вдоль оси y;
f='0'; % f - функция в правой части уравнения Пуассона;
v1=1; % v1 - параметр, значение которого определяет 
      % тип граничного условия на первой границе 
      % области х = х(1) (1 - ГУ Дирихле, 2 - ГУ Неймана); 
g1='y'; % g1 - функция в правой части граничного условия 
        % на первой границе; 
v2=1;    % на второй границе
g2='y'; 
v3=1;    % на третьей границе
g3='x'; 
v4=1;    % на четвертой границе
g4='x+1'; 
 
% Задание равномерной координатной сетки 
 
x=x0:(xn-x0)/(n-1):xn; dx=x(2)-x(1); 
y=y0:(ym-y0)/(m-1):ym; dy=y(2)-y(1); 
 
% Вычисление значений функций, заданных символьно, 
% в узлах координатной сетки 
 
F=inline(f,'x','y'); 
G1=inline(g1,'y'); 
G2=inline(g2,'y'); 
G3=inline(g3,'x'); 
G4=inline(g4,'x'); 
 
% Определение размерности СЛАУ 
 
N=n*m; 
 
% Задание матрицы коэффициентов СЛАУ размерности N x N, 
% все элементы которой равны 0 
 
a=zeros(N,N); 
 
% Задание матрицы-строки свободных членов СЛАУ 
% размерности 1 x N, все элементы которой равны 0 
 
b=zeros(1,N); 
 
% Определение коэффициентов и свободных членов СЛАУ, 
% соответствующих граничным условиям и проверка 
% корректности значений параметров v1, v2, v3, v4 
 
for j=1:m 
    b(j)=G1(y(j)); 
    if v1==1 
        a(j,j)=1; 
    elseif v1==2 
        a(j,j)=-1/dx; 
        a(j,m+j)=1/dx; 
    end
    b(m*(n-1)+j)=G2(y(j));
    if v2==1 
        a(m*(n-1)+j,m*(n-1)+j)=1; 
    elseif v2==2 
        a(m*(n-1)+j,m*(n-1)+j)=1/dx; 
        a(m*(n-1)+j,m*(n-2)+j)=-1/dx; 
    end
end
for i=2:n-1 
    b(m*(i-1)+1)=G3(x(i)); 
    if v3==1 
        a(m*(i-1)+1,m*(i-1)+1)=1; 
    elseif v3==2 
        a(m*(i-1)+1,m*(i-1)+1)=-1/dy; 
    end
    b(m*(i-1)+m)=G4(x(i));
    if v4==1 
        a(m*(i-1)+m,m*(i-1)+m)=1; 
    elseif v4==2 
        a(m*(i-1)+m,m*(i-1)+m)=1/dy; 
        a(m*(i-1)+m,m*(i-1)+m-1)=-1/dy;
    end 
end
 
% Определение коэффициентов и свободных членов СЛАУ, 
% соответствующих внутренним точкам области 
 
for i=2:n-1 
    for j=2:m-1 
     a(m*(i-1)+j,m*(i-1)+j)=-2/dx^2-2/dy^2; 
     a(m*(i-1)+j,m*(i)+j)=1/dx^2; 
     a(m*(i-1)+j,m*(i-2)+j)=1/dx^2; 
     a(m*(i-1)+j,m*(i-1)+j+1)=1/dy^2; 
     a(m*(i-1)+j,m*(i-1)+j-1)=1/dy^2; 
     b(m*(i-1)+j)=F(x(i),y(j)); 
    end 
end
 
% Решение СЛАУ
 
u=b/a'; 
 
% Преобразование вектора-строки значений искомой функции 
% в узлах координатной сетки в матрицу размерности n x m, 
% удобную для представления результатов 
% в графическом виде 
 
for i=1:n 
    for j=1:m 
        U(i,j)=u(m*(i-1)+j); 
    end
end
 
% Построение графика искомой функции U(x,y) 
 
surf(y,x,U) 
xlabel('y') 
ylabel('x') 
zlabel('U') 
grid on

График строится, все работает. Условия g1,g2,g3,g4 вроде задал логично (так мне казалось).
Препод сказал, что эти разрывы на схеме подключения означают разрывы и на графике, и что на участке V1 график должен быть в нуле.
Есть подозрение, что условия g1,g2,g3,g4 надо задавать как-то по-другому, а на нижней границе (то бишь, g3) вообще по частям. Пробовал через if, не особо что вышло.

Ну и вопрос, собственно. Как сделать эти разрывы, и возможно ли это вообще? Или я в корне ошибаюсь и решать нужно вообще иначе?

P.S.: Упомянутая выше книжка. На всякий случай
ryndin_mathphys_methods.pdf

@tynec02 · 29.05.2014, 19:22 **[ТС]**

Сделал. Правильно вроде. Мб нужно будет кому.

Matlab M

% u_xx + u_yy = 0, x,y [0,1], V1=xy, V2=x+y
 
clc
clear
close
 
x0=0; % x0 - начальная координата области решения по оси х;
xn=1; % xn - конечная координата области решения по оси х;
n=50; % n - число точек координатной сетки вдоль оси х; 
y0=0; % y0 - начальная координата области решения по оси y; 
ym=1; % ym - конечная координата области решения по оси y; 
m=50; % m - число точек координатной сетки вдоль оси y;
f=0;  % f - функция в правой части уравнения Пуассона;
 
% Задание равномерной координатной сетки  
x=x0:(xn-x0)/(n-1):xn; dx=x(2)-x(1); 
y=y0:(ym-y0)/(m-1):ym; dy=y(2)-y(1);
 
g1(1:n/2)=x(1:n/2);     % g1 - функция в правой части граничного условия 
g1(n/2:n)=x(n/2:n)*0; % на первой границе;
g2=x+1;  % на второй границе;
g3=y;    % на третьей границе;
g4=y;    % на четвертой границе;
 
% Вычисление значений функций, заданных символьно, 
% в узлах координатной сетки 
 
% Определение размерности СЛАУ
N=n*m; 
 
% Задание матрицы коэффициентов СЛАУ размерности N x N, 
% все элементы которой равны 0
a=zeros(N,N); 
 
% Задание матрицы-строки свободных членов СЛАУ 
% размерности 1 x N, все элементы которой равны 0
b=zeros(1,N); 
 
% Определение коэффициентов и свободных членов СЛАУ, 
% соответствующих граничным условиям
for j=1:m 
    b(j)=g1(j);
    a(j,j)=1;
    b(m*(n-1)+j)=g2(j);
    a(m*(n-1)+j,m*(n-1)+j)=1;
end
for i=2:n-1 
    b(m*(i-1)+1)=g3(i);
    a(m*(i-1)+1,m*(i-1)+1)=1; 
    b(m*(i-1)+m)=g4(i);
    a(m*(i-1)+m,m*(i-1)+m)=1;
end
 
% Определение коэффициентов и свободных членов СЛАУ, 
% соответствующих внутренним точкам области 
for i=2:n-1 
    for j=2:m-1 
     a(m*(i-1)+j,m*(i-1)+j)=-2/dx^2-2/dy^2; 
     a(m*(i-1)+j,m*(i)+j)=1/dx^2; 
     a(m*(i-1)+j,m*(i-2)+j)=1/dx^2; 
     a(m*(i-1)+j,m*(i-1)+j+1)=1/dy^2; 
     a(m*(i-1)+j,m*(i-1)+j-1)=1/dy^2; 
     b(m*(i-1)+j)=f; 
    end 
end
 
% Решение СЛАУ
u=b/a'; 
 
% Преобразование вектора-строки значений искомой функции 
% в узлах координатной сетки в матрицу размерности n x m, 
% удобную для представления результатов в графическом виде 
for i=1:n 
    for j=1:m 
        U(i,j)=u(m*(i-1)+j); 
    end
end
 
% Построение графика искомой функции U(x,y)
surf(y,x,U) 
xlabel('x') 
ylabel('y') 
zlabel('U') 
grid on

@Dimas_rdn · 29.05.2014, 19:53

Редкий момент, когда кто-то пишет собственное решение на свой же вопрос.

@Dimas_rdn 41 / 41 / 9 Регистрация: 22.10.2012 Сообщений: 91
	29.05.2014, 19:53	3
	Редкий момент, когда кто-то пишет собственное решение на свой же вопрос. 0