Материал: Спецглавы высшей математики. численные методы. Пантелеев И.Н

2. Выполнить следующую последовательность команд:

>>Nx=30; % число узлов по координате х

>>Nt=100; % число шагов по времени % задание начальных условий

>>j=1:Nx;

>>u(1,j)=sin(pi*(j-1)/(Nx-1))+(j-1)/(Nx-1);

% задание граничных условий

>>alpha=0;

>>beta=1;

>>i=1:Nt;

>>u(i,1)=alpha;

>>u(i,Nx)=beta;

% создание матрицы системы линейных уравнений (IX.22)

>>r=5;

>>for k=1:Nx-2 A(k,k)=1+2*r;

end;

>> for m=2:Nx-2 A(m-1,m)=-r; A(m,m-1)=-r;

end;

>>z=U1(u,A,r,Nt,Nx); % вычисление решения системы

%смешанной краевой задачи

>>surf(z); view(160,70); colormap gray

% визуализации решения смешанной краевой задачи (рис. IX.7)

Завершая обсуждение разностных методов решения дифференциальных уравнений в частных производных, отметим, что различные методы и средства для решения данного типа уравнений пакета МАТLАВ объединены в специализированный пакет Partial Differential Equations Toolboox (РDЕ Тoolbоох), подробное обсуждение которого выходит за рамки нашего курса. Заинтересованный читатель может изучить приемы работы с РDЕ Тооlbоох, ознакомившись с руководством пользователя, которое входит в состав документации, поставляемой в электронной форме вместе с пакетом МАТLАВ (файл РDЕ.pdf).

Рис. IX.7. Решение смешанной краевой задачи УЧП параболического типа

9.5. Решение уравнений с частными производными методом Монте-Карло

Рассмотрим решение задачи Дирихле уравнения Лапласа методом случайных блужданий:

uxx u yy 0, 0 < х < 1, 0 < у < 1.

Для иллюстрации метода Монте-Карло рассмотрим игру, которая называется "Блуждающий пьяница".

Правила игры:

1. Блуждания пьяницы начинаются из произвольной точки сетки (в нашем случае это точка А (рис. IX.8)).

219

Рис. IX.8. Координатная сетка, используемая в методе Монте-Карло

2.На каждом шаге пьяница случайным образом перемещается в одну из четырех соседних точек сетки (в нашем случае — одна из точек С, В, D, Е). Вероятность попадания в каждую из точек равна.

3.После перехода в соседнюю точку процесс блуждания возобновляется. Перемещения пьяного продолжаются до тех

пор, пока он не достигнет одной из граничных точек Pi. Номер граничной точки фиксируется, и на этом случайная прогулка заканчивается.

4.Производим повторение шагов 1-3 достаточное количество раз и определяем количество посещений пьяным каждой граничной точки. Отношение числа посещений пьяным каждой точки границы к полному числу испытаний определяет вероятность попадания пьяного в данную точку границы.

5.Предположим, что пьяница получает вознаграждение gi, если он достигает точки pi , и предположим, что цель игры — вычислить среднее вознаграждение R(А) для всех случайных прогулок, начинающихся из точки А . Тогда искомый средний выигрыш определяется формулой:

R( A) g1 PA ( p1 ) g2 PA ( p2 ) g12 PA ( p12 ). (IX.23)

Для чего надо играть в "Блуждающего пьяницу"? Оказывается, что среднее вознаграждение (IX.23) является решением задачи Дирихле в точке А . Данный вывод основан на двух фактах.

1. Предположим, что пьяница начал свое движение из точки, лежащей на границе. Каждая такая прогулка заканчивается в той же точке, и пьяница немедленно получает вознаграждение gi . Таким образом, среднее вознаграждение для каждой точки равно gi . 2. Теперь предположим, что прогулка начинается из внутренней точки. Тогда ясно, что среднее вознаграждение для точки К(А) будет средним арифметическим от средних вознаграждений для четырех соседних точек:

Итак, мы видим, что R(А) удовлетворяет уравнению (IX.24) в каждой внутренней точке и равно gi в граничной точке. Если gi — это значения функции g(х,у) из граничного условия в граничных точках pi, то два наших уравнения точно совпадают с двумя уравнениями, которые получены при решении задачи Дирихле методом конечных разностей. То есть величина R(А) соответствует величине uij в разностных уравнениях

где gij – значение решения в граничной точке (i, j).

Таким образом, величина R(А) действительно аппроксимирует решение уравнения с частными производными в точке А . Пример IX.5. Решение методом случайных блужданий в пакете МАТLАВ уравнения Лапласа:

1. Создайте файл Моп1е_Каг1о.т (листинг IX.3), содержащий описание функции, возвращающей численное решение уравнения Лапласа.

Листинг IX.3. Файл Monte_Carlo.m

function z=Monte_Carlo(Nx,Ny,NTrial,XL,XR,Yu,Yd) for j=1:Nx

U(1,j)=Yd(j);

U(Ny,j)=Yu(j);

end;

for i=1:Ny U(i,1)=XL(i); U(i,Ny)=XR(i);

end;

for i=2:Nx-1 for j=2:Ny-2

Xs=j;

Ys=I;

for m=1:Nx VL(m)=0; VR(m)=0;

end;

for m=1:Ny Gd(m)=0; Gu(m)=0;

end;

for k=1:NTrial x=Xs;

y=Ys;

while not(x==1)&not(x==Nx)&not(y==1)&not(y==Ny) R=rand(1,1);

if R<0.25 x=x-1;

end;

if (0.25<=R)&(R<0.5) x=x+1;

end;

if (0.5<=R)&(R<0.75)

y=y-1;

222

end;

if (0.75<=R)&(R<=1) y=y+1;

end;

end;

if x==1 VL(y)=VL(y)+1;

end;

if x==Nx VR(y)=VR(y)+1;

end;

if y==1 Gd(x)=Gd(x)+1;

end;

if y==Ny Gu(x)=Gu(x)+1;

end;

end;

s1=0;

for m=2:Nx-1 s1=s1+Gd(m)*U(1,m)+Gu(m)*U(Ny,m);

end;

s2=0;

for m=2:Ny s2=s2+VL(m)*U(m,1)+VR(m)*U(m,Nx);

end;

U(I,j)=(s1+s2)/NTrial;

end;

end;

z=U;

223

Рис. IX.9. Карта линий уровня решения уравнения Лапласа, полученного методом Монте-Карло

2. Далее необходимо выполнить следующую последовательность команд:

% задание числа узлов координатной сетки

>>Nx=14;

>>Ny=14;

% задание граничных условий

>>i=1:Nx;

>>j=1:Ny;

>>XL(i)=5;

>>Yu(j)=0;

>>XR(i)=5;

>>Yd(j)=0;

% нахождение численного решения

>> U=Monte-Carlo(Nx,Ny,XL,XR,Yu,Yd); >>contour(U,17) %визуализация численного решения (рис. IX.9)

224

Рис. IX.10. Решение уравнения Лапласа, полученное методом релаксаций

Рассмотрим эллиптическую краевую задачу в квадрате u x x (sin x)u x x 0, 0 x , 0 x ,

u(x, y) g(x, y)

на его границе.

Для того чтобы решить эту задачу, заменим uxx , uyy и sinx следующим образом:

u x x ui j 1 2ui j ui j 1, u y y ui 1 j 2ui j ui 1 j ,

sin x sin x j .

Подставим эту замену в уравнение с частными производными. Разрешив получившееся уравнение относительно uij, получаем