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

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

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

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

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

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

Matlab M    
        
            
                
                
                
                
            
        
    Скопировано
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
% 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    
        
            
                
                
                
                
            
        
    Скопировано
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
% 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

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

Новые блоги и статьи Все статьи Все блоги /
Работа с объемным DOM в javascript Htext 04.04.2025 Сегодня прочитал статью тут о расходах памяти в JS, ее утечках и т. п. И вот что вспомнил из своей недавней практики. Может, кому пригодится. Хотя, в той статье об этом тоже есть. Дело в том, что я. . .	Оптимизация производительности Node.js с помощью кластеризации run.dev 04.04.2025 Масштабирование приложений для обработки тысяч и миллионов запросов — обыденная задача для многих команд. Node. js, благодаря своей асинхронной событийно-ориентированной архитектуре, стал популярной. . .	Управление зависимостями в Python с Poetry py-thonny 04.04.2025 Стандартный инструмент для установки пакетов в Python - pip - прекрасно справляется с базовыми сценариями: установил пакет командой pip install и используешь его. Но что произойдёт, когда разные. . .	Мониторинг с Prometheus в PHP Jason-Webb 04.04.2025 Prometheus выделяется среди других систем мониторинга своим подходом к сбору и хранению метрик. В отличие от New Relic, который использует агентный подход и отправляет данные во внешнее хранилище,. . .	Пакет Context в Golang: Управление потоками и ресурсами golander 04.04.2025 Работа с горутинами в Go часто напоминает управление непослушными детьми - они разбегаются кто куда, делают что хотят и не всегда завершаются вовремя. К счастью, в Go 1. 7 появился пакет context,. . .
Контейнеризация React приложений с Docker Reangularity 03.04.2025 Контейнеризация позволяет упаковать приложение со всеми его зависимостями в автономный контейнер, который можно запустить на любой платформе с установленным Docker. Это существенно упрощает процессы. . .	Свой попап в SwiftUI mobDevWorks 03.04.2025 SwiftUI, как декларативный фреймворк от Apple, предоставляет множество инструментов для создания пользовательских интерфейсов. В нашем распоряжении есть такие API как alerts, popovers, action sheets. . .	Антипаттерны микросервисной архитектуры ArchitectMsa 03.04.2025 Хорошо спроектированная микросервисная система может выдержать испытание временем, оставаясь гибкой, масштабируемой и устойчивой к большинству проблем. Такая архитектура обладает высоким уровнем. . .	std::mutex в C++: Советы и примеры использования bytestream 03.04.2025 std::mutex - это механизм взаимного исключения, который гарантирует, что критический участок кода выполняется только одним потоком в каждый момент времени. Это простое, но могущественное средство. . .	Не удержался от оценки концепции двигателя Стирлинга. Hrethgir 03.04.2025 Сколько не пытался - она выдавала правильные схемы, причём случайно рисовала горячие области в середине, холодные по краям, трубки с краёв в низ и магнит в соединяющей, но при этой выдавала описание. . .

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