Голичев И.И. Лабораторный практикум по курсу «Численные методы»

Лабораторная работа № 11.

Решение задачи Неймана для уравнения Пуассона методом

скорейшего спуска

1. Постановка задачи

Явным методом Чебышева требуется приближенно решить однородную

задачу Неймана для уравнения Пуассона в квадрате



, которая состоит в

следующем: найти функцию

 

21

, xxU

, удовлетворяющую уравнению Пуассона

и краевым условиям

 

2 2

1 2

2 2

1 2

,f x x

x x

 

 

 

 

, (1)

0





n

U

, (2)

где



– граница квадрата



,

n

– внешняя нормаль к



.

Функция, удовлетворяющая краевым условиям (2)

32

3

1

3

2

1

xxxx

U









.

Вычислим

 

21

2

1

2

2 xx

x

U

x

U









.

Возьмём по определению в качестве правой части уравнения (1)

   

1 2 1 2

, 2 .f x x x x 

.

2. Теоретическая часть

Для аппроксимации производной в граничных условиях (2) используем

формулы второго порядка точности через центральные разности. Для этого

вводим дополнительные точки за пределами сеточной области

       

Njihihxhihxjhhxjhhx

NiijNjh

,...,1,0,,1,,,,,1,,:

1,1,,1,1





.

Получаем:

,0

2

,0

2

1,1,1,1,











h

yy

h

yy

NiNiii

(3)

,0

2

,0

2

,1,1,1,1











h

yy

h

yy

jNjNjj

Nji ,...,1,0, 

.

На расширенной сетке аппроксимируем задачу (1), (2) в виде

1, , 1, , 1 , , 1

,

2 2

i j i j i j i j i j i j

h i j

y y y y y y

f

h h

   

   

   

,

Nji ,...,1,0, 

. Пользуясь соотношениями (3), дополнительные неизвестные

1,1,,1,1

,,,

 NiijNj

yyyy

исключаем. Тогда получим

31

ji

jijijijijiji

h

f

h

yyy

h

yyy

,

2

1,,1,

2

,1,,1

22













,

j

jjjijj

f

h

yyy

h

yy

,0

2

1,0,01,

2

,1,0

1

222













,

jN

jNjNjNjNjN

f

h

yyy

h

yy

,

2

1,,1,

2

,1,

2

222













,

0,

2

0,10,0,1

2

1,0,

3

222

i

iiiii

f

h

yyy

h

yy













,

(4)

Ni

NiNiNiNiNi

f

h

yyy

h

yy

,

2

,1,,1

2

1,,

4

222













,

0,0 1,0 0,1 ,0 1,0 ,1

0,0 0,0 ,0 ,0

2 2

4 2 2 4 2 2

,

N N N

N N

y y y y y y

f f

h h



   

     

,

0, 0, 1 1, , 1, , 1

0, 0, , ,

2 2

4 2 2 4 2 2

,

N N N N N N N N N

N N N N N N

y y y y y y

f f

h h

  

   

     

,

, 1,..., 1.i j N 

Обозначим через



пространство функций

y

, заданных на сетке

h



со

скалярным произведением

 

2

, ,

0 0

, .

N N

i j i j

i j

y v h y v

 





Сеточную функцию

v

будем рассматривать как вектор с координатами

 

Njiv

ij

,...,1,0, 

. В пространстве



определим оператор



, который сеточной

функции

y

с координатами

 

Njiy

ij

,...,1,0, 

сопоставляет сеточную функцию



с координатами

ij



, где

ij



определяется левыми частями уравнений (4). То

есть

1

2

3

4

0,0

,0

0,

,

, 1 N-1, 1 j N-1,

, 0, 1 1,

, , 1 1,

, , 1 1, 0,

( ) , , 1 1, ,

, 0, 0,

, , 0,

, 0, ,

, , .

h

ij ij

N

N N

если i

если i j N

если i N j N

если i N j

Ay w если i N j N

если i j

если i N j

если i j N

если i N j N



    



    



    



    



      





  



  



  



  



Таким образом, краевая задача (1), (2) аппроксимируется операторным

уравнением

fy 

, (5)

где

f

– сеточная функция,

 

ijij

fxf 

. Заметим, что оператор



симметрический

и однородное уравнение

32

0y

(6)

имеет нетривиальные решения

consty 

на сетке

h



, кроме того, любое

решение (6) есть

const

на

h



. Обозначим





подпространство сеточных

функций из



, ортогональных

const

(например ортогональных к

1

0

y

на

h



). Такую функцию будем обозначать



1

.

V

принадлежит





тогда и только

тогда, когда

 









N

ji

vhv

0,

,

2

01,0

.

Заметим, что

y

принадлежит



при любом

y

из



.

Действительно, поскольку

 

01 



, то

    





yyy ,01,1,0

.

Отсюда следует, что уравнение (5) имеет решение тогда и только тогда,

когда

f

принадлежит





и любое решение определено с точностью до

const

на

h



.

Условимся выбирать решение уравнения (5) принадлежащее





. Если

y

какое-то решение (5), то



 1



y

принадлежит





, где

 

2

0

1







N

y

N

ij



,

а



1



– сеточная функция, равная



на всех точках сетки. Действительно

   

2 2

1 ,1 1 1N

  



  

, поэтому

   

 

2

, 0

2

, 1

,1 1 ,1 1 0

1

N

ij

N

i j

ij

i j

y

y y N

N





  





    





.

Уравнение (5) на подпространстве





невырожденное и

 



 yyyyy ,,

2

max

2

min



,

где

max



– наибольшее собственное значение оператора



,

min



– наименьшее

неравное нулю собственное значение. Если для решения уравнения (5) мы

воспользуемся итерационным процессом, то нужно следить, чтобы итерации

не выходили из





. Этого легко добиться при использовании явного

итерационного процесса

1

, 0,1...

k k

k

y y

y f k

T





   

В этом случае

1k k k

k k

y y T y T f



   

(7)

Откуда следует, что если

0

, , то

k

y f y

  

  

.

Заметим здесь, что если для решения задачи (5) потребуется большое

число итераций, то, в силу накопления погрешностей

k

y

может выходить из





, поэтому следует через некоторое количество итераций подправлять

k

y

, то

есть заменять

k

y

на

33

 

, 0

2

1 ,

1

N

k

ij

i j

k k

y

y y

N

 





  







. (8)

3. Алгоритм

1) Задаём

10, 0,1N h 

, начальное приближение

 

0 0

0 , 0, .

ij

y y i j N  

2) Вычисляем параметры

k

T

по формуле

( , )

, .

( , )

k

k k

r r

T r Ay f

Ar r

  

3) Вычисляем

1k

y



по формуле (7).

4) Через 10 итераций подправляем

k

y

по формуле (8).

5) Вычисление проверить до совпадения первых четырех знаков в последних

итерациях.

4. Оформление результатов работы.

Результаты вычислений представить в виде трех таблиц: две

последовательные итерации с совпадением первых четырех знаков и значение

точного решения на сетке.

34

Лабораторная работа № 12.

Решение задачи Дирихле для уравнения Пуассона

методом переменных направлений

1. Постановка задачи

1) Методом переменных направлений решить задачу Дирихле для уравнения

Пуассона

 

21

2

1

2

, xxf

x

U

x

U









(1)

0



U

, (2)

в области

  

10,10:,

2121

 xxxxD

, где



– граница квадрата

D

. Пусть

 

2

11221

2, xxxxxxf 

, (3)

тогда

   

2211

11 xxxxU 

– точное решение задачи (1), (2).

Написать программы для реализации метода переменных направлений на

любом языке программирования. Оценить точность решения задачи.

2. Задача Дирихле для уравнения Пуассона

Пусть

  

10,10:,

2121

 xxxxD

– замкнутый квадрат,



– его

граница,

 

21

, xxf

– заданная на

D

достаточно гладкая функция. Задача

Дирихле состоит в следующем: требуется найти непрерывную на

D

функцию

 

21

, xxU

, удовлетворяющую на открытом квадрате

D

уравнению Пуассона

(1) и обращается в 0 на границе квадрата.

Задача Дирихле (1), (2) имеет единственное решение

 

21

, xxU

. Положим

Nh 1

,

 

mkkmmk

xxffmhxkhx

2121

,,, 

. Построим сетки

  

Nmkxxw

mkh

...1,0,:,

21



,

  

1...2,1,:,'

21

 Nmkxxw

mkh

,

hhh

www '\

*



(

*

h

w

– множество узлов, лежащих на



).

Заменим исходную дифференциальную задачу (1), (2) разностной задачей

1, , 1, , 1 , , 1

2 2

, , 1, 2... 1,

k m k m k m k m k m k m

km

V V V V V V

f k m N

h h

   

    

(4)

0

,



mk

V

на

*

h

w

. (5)

Введём обозначения:

2

,1,,1

,1

2

h

VVV

V

v

mk

v

mk

v

mk

v

mk





,

(6)

2

1,,1,

,2

2

h

VVV

V

v

mk

v

mk

v

mk

v

mk





,

(7)

v

mk

v

mk

v

mk

VVV

,2,1,



.

35

Таким образом, наше уравнение (4) можно переписать в виде

mk

v

mk

v

mk

fVV

,,2,1



1...2,1,  Nmk

.

(4*)

3. Метод переменных направлений для уравнения Пуассона

Сопоставим задачу (1), (2) с родственной её нестационарной задачей о

распространении тепла.

Пусть

 

 

 

 

 

 

1 2 1 2

, , : 0 1, 0 1, 0 ,

, : 0 1, 0 1 ,

' , , : 0 1, 0 1, 0 ,

'

R x x t x x t T

D x x x x

R x x t x x t T

S R R

      

    

      



То есть

S

– множество граничных точек прямоугольного

параллелепипеда

R

, не принадлежащих

'R

,

 

txxf ,,

21

– заданная на

R

достаточно гладкая функция.

Требуется найти непрерывную на

R

функцию

 

txxW ,,

21

,

удовлетворяющую на

'R

уравнению теплопроводности

 

21

2

1

2

, xxf

x

W

x

W

t

W















 

 Dyx,

,

0t

(8)

и кроме того, подчиняющуюся на

S

дополнительному условию

0W

на

S

. (9)

Условие (9) включает в себя как начальное условие

 

00,,

21

xxW

при

0t

, так и однородные краевые условия первого рода.

Смешанная задача (8), (9) имеет единственное решение

 

txxW ,,

21

.

В задаче (8), (9) источник тепла

 

21

, xxf

и температура на границе

S

не

зависят от времени

t

. Естественно ожидать поэтому, что при

t

будет

выполнятся соотношение

0



t

W

, откуда

   

2121

,,,lim xxUtxxW

t





,

поэтому можно предположить, что для достаточно больших значений

t

,

например для

Tt 

, будет с необходимой точностью верно приближённое

равенство. Здесь положена в основу идея о стабилизации при

t

решения

уравнения теплопроводности к решению уравнения Пуассона, если

f

не

зависит от

t

, т.е.

   

2121

,,, xxUtxxW 

.

На этой закономерности основана идея метода решения стационарных

задач, состоящая в замене их подходящими нестационарными задачами. Этот

метод и ряд его модификаций принято называть методом установления.

Запишем для задачи (8), (9) простейшую двухслойную разностную схему

1

, ,

, , 1,2,... 1, 1,2,...,

v v

k m k m

v

k m k m

V V

V f k m N v





     

0

,



v

mk

V

на



36

и двухслойную разностную схему переменных направлений или дробных

шагов

mk

v

mk

v

mk

v

mk

v

mk

fVV

VV

,

1

,2

21

,1

1

,

21

,

2









, (10)

mk

v

mk

v

mk

v

mk

v

mk

fVV

VV

,,2

21

,1

21

,,

2







, (11)

,...2,1,1,...2,1,  vNmk

.

В разностной схеме (10), (11) шаг



по времени делится на два полушага.

Разностное уравнение (10) отвечает первому полушагу, в нём величины

1

,

v

mk

V

и

1

,2





v

mk

V

считаются уже известными (в частности,

NmkV

mk

,...1,0,,0

0

,



), а

неизвестные имеют верхний индекс

21v

. Правая часть задана. Перепишем

разностное уравнение (10), предварительно умножив его на

2





, следующим

образом:

21

,

21

,1

2

21

,

2

21

,1

2

1

2

























v

mk

v

mk

v

mk

v

mk

FV

h

V

h

V

h



, (12)

где

 

1

,,

1

,2

21

,

2





v

mkmk

v

mk

v

mk

VfVF



известно, и присоединим к разностному уравнению краевые условия

0

21

,0



v

m

V

,

0

21

,



v

mN

V

(13)

в соответствии с условием (9).

Разностная задача (12), (13) распадается на

1N

независимых

трёхточечных разностных краевых задач, отвечающих каждому

фиксированному значению

m

,

1,...2,1  Nm

. Разностная краевая задача (12),

(13) решается методом прогонки при каждом

m

отдельно. Прогонка

осуществляется по индексу

k

, то есть в направлении оси

1

x

.

После того как найдены все неизвестные

21

,

v

mk

V

на промежуточном слое с

номером

21v

, переносим их в разностном уравнении (11), соответствующем

второму полушагу, вправо. Это разностное уравнение переписываем в виде

v

mk

v

mk

v

mk

v

mk

FV

h

V

h

V

h

,1,

2

,

2

1,

2

1

2



















, (14)

где

 

21

,,

21

,1,

2





v

mkmk

v

mk

v

mk

VfVF



известно, и присоединяем к уравнению (14) в соответствии с условием (9)

краевые условия

0

0,



v

k

V

,

0

,



v

Nk

V

. (15)

Задача (14), (15) тоже не распадается на

1N

независимых трёхточечных

разностных краевых задач, отвечающих каждому фиксированному

1,...2,1,  Nkk

. Каждая такая задача решается методом прогонки. Прогонка

осуществляется теперь уже по индексу

m

, то есть в направлении оси

2

x

.

4. Метод прогонки

37

Рассмотрим систему линейных алгебраических уравнений следующего

специального вида:

jjjjjjj

FzBzCzA 

 11

,

,1,...2,1  Mj

MMMM

vzazvzaz 

 10100

,

где

M

zzz ,...,

10

– неизвестные,

0 0

, , , , , , ,

j j j j M M

A B C F a v a



– заданные числа,

причём

1,1,0

0



Mjjjj

aaABAC

.

Решение ищем по формуле

kkkk

vzaz 

1

, (*)

где

1

11

1

,









kkk

kk

k

kkk

k

aAC

FvA

v

aAC

B

a

. (**)

Таким образом, коэффициенты уравнений (*), связывающих соседние

значения

1,...2,1,,

1





Mkzz

kk

можно определить из рекуррентных

соотношений (**), поскольку

0

a

,

0

v

заданы.

Находим неизвестную

M

z

по формуле

1







MM

MMM

M

aa

vav

z

(***)

и в обратном порядке находим все неизвестные

021

,..., zzz

MM 

. Процесс

вычисления коэффициентов

1,2,1,,  Mkva

kk

по формулам (**) называется

прямым ходом прогонки, а нахождение неизвестных

, , 1,...0

m

z m M M 

, по

формулам (***), (*) – обратным ходом прогонки.

5. Алгоритм решения задачи Дирихле для уравнения Пуассона

Как уже было замечено выше, будем искать решение задачи Дирихле для

уравнения Пуассона в области

  

10,10:,

2121

 xxxxD

1) начальные условия

N

– натуральное число,

Nh 1

шаг по

1

x

и по

 

mkkmmk

xxffmhxkhxx

21212

,,,, 

.

,...2,1,0,1,...2,1.0,0,0,0

,0,,,0

 vNkVVVV

v

Nk

v

k

v

mN

v

m

0

,

0

k m

V 

1,...2,1,  Nmk

– начальное приближение. Полагаем

2sin( )h h

 



2) Прогонка в направлении оси

1

x

Решим методом прогонки при каждом фиксированном

1...2,1  Nm

систему уравнений (12)

21

,

21

,1

2

21

,

2

21

,1

2

1

2

























v

mk

v

mk

v

mk

v

mk

FV

h

V

h

V

h



,

где

 

1

,,

1

,2

21

,

2





v

mkmk

v

mk

v

mk

VfVF



известно,

1,...2,1,

2

,1,

2

222















 Nk

h

B

h

C

h

A

kkk



.

,...3,2,1,0,0,0,0

21

,

21

,

21

,0

21

,0





vvava

v

mN

v

mN

v

m

v

m

38

Перепишем уравнение (12) в виде:

21

,

21

,1

21

,

21

,1











v

mk

v

mkk

v

mkk

v

mkk

FVCVBVA

, (12*)

где

 

1

,,

1

,2

21

,

2





v

mkmk

v

mk

v

mk

VfVF



.

При

1v

получаем

211

,

211

,1

211

,

211

,1











mkmkkmkkmkk

FVBVCVA

(12*),

где

mkmk

fF

,

211

,

2







Прогонка осуществляется при каждом фиксированном

1,...2,1  Nm

.

Прямой ход прогонки

Вычислим коэффициенты

1...2,1,

21

,

21

,1





Nkva

v

mk

v

m

21

,011

21

,

21

,01

21

,1

21

,011

1

21

,1













v

m

v

mk

v

m

v

m

v

m

v

m

aAC

FvA

v

aAC

B

a

.

Так как

0,0

21

,0

21

,0



 v

m

v

m

va

, то получаем

1

21

,

21

,1

1

21

,1

C

F

v

C

B

a

v

mk

v

m

v

m









и далее

21

,1

21

,

21

,1

21

,1

21

,

21

,













v

mkkk

v

mk

v

mkk

v

m

v

mkkk

k

v

mk

aAC

FvA

v

aAC

B

a

.

После того, как будут вычислены коэффициенты

1...2,1,

21

,

21

,





Nkva

v

mk

v

mk

, вычислим

21

,

v

mN

z

21

,1

21

,

21

,1

21

,

21

,

21

,

1











v

mN

v

mN

v

mN

v

mN

v

mN

v

mN

va

vav

V

.

Так как

0,0

21

,

21

,



 v

mN

v

mN

va

, то получаем

1 2

,

0

v

N m

V





.

При

1v

получаем

,

2

,

1

,

211

,1

1

211

,1

C

f

v

C

B

a

mk

mm







1...2,1,

2

,

211

,1

,

211

,1

211

,

211

,1

211

,















Nk

aAC

f

vA

v

aAC

B

a

mkkk

mk

mkk

mk

mkkk

k

mk



1 1 2

,

0.

N m

V





Обратный ход прогонки

После того как будут найдены все

1 2 1 2 1 2

, , ,

, , 0, 1,2... 1

v v v

k m k m N m

a v V k N

  

  

найдём

все неизвестные

1...1,

21

,





NkV

v

mk

по формуле

1,...,2,1,

21

,

21

,1

21

,

21

,









NkvVaV

v

mk

v

mk

v

mk

v

mk

.

Таким образом, вычисляются

1 2 1 2 1 2

, 0, ,

, 1,2,... 1; 0, 0

v v v

k m m N m

V k N V V

  

   

из

начальных условий.

39

При

1v

получаем

1 1 2 1 1 2 1 1 2 1 1 2

, , 1, ,

, 1,2,... 1

k m k m k m k m

V a V v k N

   



   

3) Прогонка в направлении оси

2

x

.

Решим методом прогонки при каждом фиксированном

1...2,1  Nk

систему

уравнений (14)

v

mk

v

mk

v

mk

v

mk

FV

h

V

h

V

h

,1,

2

,

2

1,

2

1

2



















, (14)

где

 

21

,,

21

,1,

2





v

mkmk

v

mk

v

mk

VfVF



известно из предыдущих вычислений,

1,...2,1,

2

,

2

1,

2

222















 Nm

h

B

h

C

h

A

mmm



,

,...3,2,1,0,0,0,0

,,0,0,

 vvava

v

Nk

v

Nk

v

k

v

k

в силу граничных условий.

Перепишем систему уравнений (14) в виде:

v

mk

v

mkm

v

mkm

v

mkm

FVBVCVA

,1,,1,





,

(14*)

где

 

21

,,

21

,1,

2





v

mkmk

v

mk

v

mk

VfVF



,

k

v

Nk

v

Nk

v

Nk

v

k

v

k

v

k

vVaVvVaV 

 1,,0,1,0,0,

,

При

1v

получаем

1 1 1 1

, 1 , , 1 ,m k m m k m m k m k m

A V C V B V F

 

  

,

где

 

1 1 1 2 1 1 2

, 1 , , ,

2

k m k m k m k m

F V f V



 

   

Прогонка осуществляется при каждом фиксированном

1,...2,1  Nk

.

Прямой ход прогонки

Вычисляем коэффициенты

1...2,1,,

,,

 Nmva

v

mk

v

mk

v

k

v

mk

v

k

v

k

v

k

v

k

aAC

FvA

v

aAC

B

a

0,11

,0,1

1,

0,11

1

1,









Так как

0,0

0,0,



v

k

v

k

va

, то получаем

1

,

1,

1

1,

C

F

v

C

B

a

v

mk

v

k

v

k





v

mkkm

v

mk

v

mkm

v

mk

v

mkmm

m

v

mk

aAC

FvA

v

aAC

B

a

1,

,1,

,

1,

,









После этого, как будут вычислены коэффициенты

1...2,1,,

,,

 Nmva

v

mk

v

mk

вычислим

v

Nk

z

,

v

Nk

v

Nk

v

Nk

v

Nk

v

Nk

v

mN

va

vav

V

1,,

1,,,

,

1







.

40

Голичев И.И. Лабораторный практикум по курсу «Численные методы»

Подождите немного. Документ загружается.