УП_Вычисл_матем_Кузина-Кошев

F tn 1, u tn 1	dF		o	2	.
	F tn , u tn
		dt tn

Получим:

1

du

2

d 2u

o

3

2

dt

tn

dt 2

tn

1

F tn

, u tn

1

F tn , u tn

dF

o

2

2

dt

2

2

tn

du

F

tn , u tn o

2

,

dt tn

так как

d 2u

dF

. Следовательно,

точное решение удовлетворяет разно-

dt2

dt

стному уравнению с точностью до o 2 .

Метод решения уравнения (110), использующий схему (117), не позволяет явно выразить un 1 через un . Следовательно, для решения

системы уравнений (117) необходимо применять итерационный процесс, эффективность которого основывается на хорошем начальном приближении к решению.

На практике такой же порядок точности, как и в методе (117), можно получить двумя итерациями следующего вида:

1)по формуле Эйлера вычислить u~n 1 un F tn ,un ;

2)подставить u~n 1 в правую часть схемы (117):

					~
	un 1 un				F tn ,un F tn 1, un 1	.
				2
	Следовательно,
	un 1 un	1	F tn ,un F tn 1, un F un ,tn .
		2

Такая разностная схема дает явное выражение для un 1 через un .

Примечания:

~

1. Разность (un 1 un 1 ) позволяет судить о величине шага и служит осно-

ванием для автоматического выбора шага.

2. Обобщение такого способа построения разностных уравнений позволяет получить разностные методы более высоких порядков точности.

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

~

межуточных значений u , вычисляемых последовательно, получим методы типа

Рунге – Кутта.

141

11.3.Схемы Рунге – Кутта

1)Схема второго порядка.

	На первом этапе находим промежуточное значение
			u	n un F tn ,un .
	На втором этапе определяем значение un 1 по формуле
	un 1 un 1 F tn ,un F tn ,				u	n .	1
	Такая схема имеет порядок аппроксимации o 2 при
							2 .
	В частном случае,	при 1 , 1, получаем схему, называемую
		2
	«предиктор – корректор» (счет – пересчет):
	un 1 un	1 F tn ,un F tn 1,un F tn ,un .
		2

Заметим, что такая схема уже была рассмотрена (см. подразд. 11.2.2): сначала схема Эйлера с шагом (предиктор), затем схема с полусуммой (корректор).

Идея метода «предиктор – корректор» часто используется при написании разностных схем для уравнений математической физики с частными производными.

2) Схема четвертого порядка. Она имеет вид

un 1 un 1

k1 un

2k2 un 2k3 un k4 un ,

n 0,1, ; u 0 u0 ,

6

k

F t

;

где k

n

, u

n

k

F t

n

,

u

n

1

;

2

1

2

2

k

2

F t

.

k

3

F t

n

,

u

n

; k

4

n

,

u

n

k

3

2

2

При определении un 1

по заданному un необходимо четыре раза

вычислить функцию F t,u

.

Вычисления производятся последовательно.

Для невязки (точности аппроксимации) получим выражение

n

1

k1 u tn 2k2

u tn 2k3 u tn k4 u tn u tn 1 u tn

.

6

С помощью разложения в ряд Тейлора в окрестности точки tn

можно

получить оценку

n

o 4 ,

то

есть

схема имеет четвертый порядок

аппроксимации.

Предположим, чтоужевычисленыn значений u0 , u1 , , un . Тогдаможно

построить

интерполяционный

полином

Pn t

по табличной функции

142

ui , F ti ,ui , где i n, n 1, , n k 1, и при вычислении un 1 использовать
этот полином Pn t вместо функции	F t	на интервале tn , tn 1
(экстраполяция):
tn 1	t dt.
un 1 un Pn

tn

ПолучимпринциппостроенияметодовтипаАдамса. Неудобствометода заключается в том, что для его применения необходимо иметь n точек u1 , u2 , , un , построенных другим методом.

Если F t,u имеет особенность в некоторой точке, необходимо в окрестности этой точки аппроксимировать F t,u гладкой функцией и

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

На ошибку аппроксимации накладывается ошибка округления. Неточность, возникающая при отдельном шаге из-за округлений, не должна превосходить неточность аппроксимации. В этом случае вместе с

уменьшением шага интегрирования необходимо пропорционально увеличивать точность вычислений.

11.4. Разностные методы решения ОДУ второго порядка

Для получения разностного уравнения на примере линейного дифференциального уравнения второго порядка необходимо:

–заменить область непрерывного изменения аргументов дискретным множеством точек (сеткой);

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

Рассмотрим линейное дифференциальное уравнение

		u p(x)u q(x)u f (x) ,	(118)
где функции p(x), q(x),		f (x) непрерывны.
Требуется найти решение этого уравнения, удовлетворяющее краевым
условиям:
u 0 0 ,	u 1 1 ,	где 0 , 1 – заданные числа;
q x 0,	x 0,1 .
Это частный случай краевой задачи.			x xi (
Дифференциальное		уравнение (118) во внутренних точках

i 1, 2, , N 1) интервала [0, 1] приближенно можно заменить системой линейных уравнений

		ui 1 2ui ui 1		pi		ui 1 ui 1	qiui fi ,	(119)
		h2				2h
где i 1, 2, , N 1,	u0 0 ,		uN 1 ,		pi p xi ,		qi q xi ,	fi f xi .

143

Если принять

b

1

pi

, c

i

2

q

,

a

i

1

pi

,

i

h2

2h

h2

i

h2 2h

fi fi ,

1 2 0,

1 0 ,

2 1,

то получим краевую задачу для разностного уравнения 2-го порядка с переменными коэффициентами, аналогичную задаче (96):

aiui 1 ciui biui 1 fi ,

u0 1u1 1, uN 2uN 1 2 .

Эту задачу удобно решать методом прогонки (см. подразд. 10.1.1). Отметим, что порядок аппроксимации краевой задачи разностной

схемой равен o h2 и разностная схема устойчива. Пример 4 2

Решить	уравнение y 2xy 2 y 4x при краевых условиях:
	0; y(1) 1 e 3,718.
y(0) y (0)

Решение.

Решение получено в ячейках MS Excel (рис. 33).

Выбрав шаг h 0,1, заменим исходное уравнение и краевые условия системой конечно-разностных уравнений согласно формуле (119).

y

i 2

2y

i

1

y

i

2xi

y

i 1

y

i

2yi 4xi ,

i 1,2, ,8;

0,01

0,1

y

y

0

y0

1

0;

0,1

y

3,718.

10

После приведения подобных членов получим

yi 2 2 0,2xi

yi 1

0,98 0,2xi yi 0,01 4xi .

Заносим в ячейки A2–G2 исходные

данные:

a0

1; a1 1;

A 0 –

коэффициенты начального условия;

b0 1;

b1 0;

B 3,718 – коэффициен-

ты конечного условия и шаг h 0,1.

Обозначения:

y

m y

k

y

h2 f

;

i 2

i 1

i

i

i

i

a

0

y

0

a

y1 y0

A;

b y

n

b

yn yn 1

B.

1

h

0

1

h

Заполняем расчетную таблицу.

В ячейки B5–B15 заносим xi (от 0 до 1 с шагом 0,1), в ячейки A5–A15 – соответствующие значения индексов i.

144

Рис. 33. Пример решения дифуравнения методом прогонки
Прямой ход
В принятых обозначениях mi 2 0,2xi ;	ki 0,98 0,2xi ;	fi 4xi .

Заполним ячейки С5–С13 значениями mi, ячейки D5–D13 значениями ki, ячейки–E13 значениями fi.

В ячейки F5 и G5 введем соответственно формулы:

	c0				a1 a0h	; d0	k0 Ah	f0h2 .
			m0	a1	a0h k0a1		a1 a0h

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

ci mi 1ki ci 1 ; di fi h2 ki ci 1di 1.

	Обратный ход		1cn 2dn 2				Bh
	В ячейку H15 вводится формула	yn
									.
				1 c	n 2			h
			1				0

Ячейки H14–H6 заполняем, применяя рекуррентные формулы

145

yn 1 cn 2 dn 2 yn .

Ячейка H5 вычисляется по формуле y0 a1 y1 Ah . a1 a0h

Для сравнения с точным значением последний столбец вычисляется по формуле y x ex2 .

Контрольные вопросы и задания

1.В каком узле определяется производная при численном дифференцировании в случае использования трех и пяти узлов?

2.Как численно определить значение производной в произвольной

точке?

3.Какие погрешности возникают при численном дифференцировании таблично заданной функции?

4.Получите общую формулу приближенного вычисления второй

производной функции. Примените её для вычисления f (x) в точке х = 4, если f (x) 2x2 . Сравните с точным значением.

5.В каких случаях результаты численного дифференцирования по формулам при 3-х и 5 узлах, окружающих точку дифференцирования, совпадают? Покажите, что в этом случае совпадают формулы дифференцирования.

6.Что является решением дифференциального уравнения?

7.Какие методы можно применять для решения системы алгебраических уравнений, полученных при аппроксимации дифференциального уравнения второго порядка разностной схемой?

146

12. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

ВЧАСТНЫХ ПРОИЗВОДНЫХ

12.1.Уравнения математической физики

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

Необходимость решения таких уравнений возникает, когда исследуется объект или процесс, свойства которого описываются функциями нескольких переменных (например, времени и расстояния).

Дифференциальные уравнения в частных производных называются

уравнениями математической физики.

В общем виде нелинейное уравнение в частных производных может

быть записано в виде
					(n) (n)	(n)
F x, t, u, ux , ut		, uxx , utt		, utx , ..., ux...x , ut...t		, ut...x 0,

где t, x – независимые переменные; u – искомая функция; ux

ut	u(x,t)	,	uxx	2u(x,t)	,	utt	2u(x,t)	,	utx	2u(x,t)
	t			x2			t2			t x

u(x,t) ,

x

,…,

ux( n...)x	nu(x,t)	,	ut(...n t)	n u(x,t)	,	ut(...n )x	nu(x,t)	,	n k l, – частные
	xn			t n			t k xl

производные первого, второго и n-го порядков соответственно по переменным t и x.

Решением уравнения называется функция u u(t, x) , имеющая частные

производные до требуемого порядка и обращающая это уравнение в тождество. График решения представляет собой поверхность в пространстве Otxu (интегральная поверхность).

Порядком уравнения называется наивысший порядок частных производных, входящих в уравнение.

Уравнениеназываетсялинейным, еслифункцияF линейнаотносительно искомой функции u и ее производных.

Если коэффициенты уравнения не зависят от x и t, то уравнение называется уравнением с постоянными коэффициентами, в противном случае – уравнением с переменными коэффициентами.

Уравнение называется однородным, если правая часть равна нулю для всех x и t, в противном случае – неоднородным.

147

Рассмотрим вопросы сходимости, аппроксимации и устойчивости на примере уравнения в частных производных первого порядка

		U		a U		F x,t ,
		t			x
	x , t 0, U x,0 x .
Известно аналитическое решение этого уравнения:
					t			dt		.
U x,t x at F x at at ,t
				0
Построим численный метод решения задачи сведением ее к системе
конечно-разностных уравнений. Для этого введем сетку
x xk kh, k 0, 1, ;						h шагпокоординатеx ,
t tn n ,		n 0,1, ;				шагпокоординате t .
Вместо функций U x,t ,		F x,t , x будем рассматривать сеточные функ-
ции, соответственно uk ,n , fk ,n ,				k . Заменив производные в исходном урав-
нении их разностными аналогами F (xk , tk ) fk ,n , k (xk ), получим:
	uk ,n 1 uk ,n			a	uk 1,n uk ,n		fk ,n , uk ,0			k .	(120)
						h

Заметим, что такой способ аппроксимации непрерывного уравнения конечно-разностным аналогом не является единственным.

Вычислительный алгоритм заключается в последовательном нахождении решения на каждом слое по переменной t (слоем называются узлы, лежащие на одной прямой t tk при фиксированном k 0,1, , m ).

uk ,n 1 1 a h uk ,n a h uk 1,n fk ,n .

12.1.1. Аппроксимация

Прирешениизадачиразностнымметодомнадознать, скакойточностью

решение разностной задачи приближает решение исходной задачи

uk ,n U xk ,tn

uk ,n U k ,n uk ,n .

Подставив выражения для ukn

в уравнение (120), получим:

u

k ,n 1

u

k ,n

a

u

k 1,n

u

k ,n

f

u

k ,n 1

u

k ,n

a

u

k 1,n

u

k ,n

k ,n

.

(121)

h

h

Разложим в ряд Тейлорафункцию U x,t вточке xk ,tn попеременным

x и t до величин второго порядка малости:

u

u

U

ο(

2

);

k ,n 1

k ,n

t k ,n

uk 1,n uk ,n

U

ο(h

2

).

h

t

k ,n

148

Отсюда следует, что

	uk ,n 1 uk ,n				a	uk 1,n uk ,n			fk ,n
							h		(122)
		U		U
			a				F xk ,tn ο( ,h).
		t
				x k ,n

Так как	U		a	U	fk ,n , получаем точность аппрокси-
	F xk ,tn	t
				x k ,n

мации исходного дифференциального уравнения конечно-разностным аналогом, равную o ,h .

12.1.2. Сходимость

Подставив выражение (122) в равенство (121), после преобразований получим:

o ,h .

u

k ,n 1

1 a

u

k ,n

a

u

k 1,n

(123)

h

h

Выберем и h,

a , 0 a

1

. Тогда

h

o ,h

u

k ,n 1

1 a

u

k ,n

a

u

k 1,n

h

h

max

uk ,n

,

uk 1,n

o ,h .

k

o ,h .

Обозначим

un

max

uk,n

. Отсюда

un 1

un

k

Следовательно, за N шагов по τ получим:

uN

u0

N o ,h .

Значит, если t N const τ, h 0, N ,

так как u0 k k 0 ,

то uN o ,h . Следовательно, ошибка аппроксимации стремится к 0 при,h 0 , т.е. разностная схема является сходящейся.

12.1.3. Устойчивость

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

uk ,n 1 k .

Подставив ее в выражение (123), получим:

k

k 1

k

u

k ,n 1

1

a

1

a

1

1

2a

1

,

h

h

h

149

2a

. За N слоев она будет

т.е. ошибка накапливается с множителем

1

иметь значение

h

N

k

u

k ,n N

1 2a

1 .

h

2a

1

ошибка уменьшается с ростом N, в

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

1

h

противном случае она экспоненциально возрастает.

Контрольные вопросы и задания

1.Какие уравнения называют уравнениями математической физики?

2.Что такое порядок дифференциального уравнения в частных производных?

3.Какое дифференциальное уравнение в частных производных называется линейным, однородным?

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

5.Единственным ли способом можно аппроксимировать непрерывное дифференциальное уравнение в частных производных конечно-разностным аналогом?

6.Что называют слоем по переменной t?

7.В каких случаях ошибка аппроксимации дифференциального уравнения конечно-разностным аналогом возрастает, а в каких – уменьшается?

150