Выч. мат. Лекции. 6 семестр
.pdfфункция p(x) n(x)Qm(x) не меняет своего знака на отрезке [a; b] и интеграл от этой функции не равен 0. С другой стороны этот интеграл должен быть равен 0, так как
|
n(x) ортогональна на [a; b] любому полиному степени m < n. Остается принять, что |
||||||||||||||||||||||||||||||
m = n, т.е. что все нули полинома |
n(x) принадлежат [a; b], что возможно только в |
||||||||||||||||||||||||||||||
случае, когда все нули простые. |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||
|
Построение квадратурной формулы. |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||
b |
1. Строим полиномы |
n(x) â âèäå |
|
n(x) = xn + an 1xn 1 + + a1x + a0 из условий |
|||||||||||||||||||||||||||
Ra |
p(x) n(x)xj dx = 0 äëÿ j = 0; 1; : : : ; (n 1). |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||||
|
|
|
b |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Обозначим числа Ra |
p(x)xj dx = b j. Будем предполагать, что эти значения могут |
|||||||||||||||||||||||||||||
j = 0; 1; : : : ; (n 1), òî |
Ra |
|
|
|
|
|
|
n 1 |
|
|
|
|
|
|
|
|
|
|
|
1 |
0 |
|
|
|
|||||||
быть вычислены точно. Так как p(x)[xn + a |
|
xn |
|
|
1 |
+ |
|
+ a |
x + a |
]xj dx = 0 ïðè |
|||||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
a0 0 + a1 1 + + an 1 n 1 + n = 0; (j = 0); |
|
|
|
|||||||||||||||||||||||||
|
|
|
a0 1 + a1 2 + + an 1 n + n+1 = 0; (j = 1); |
|
|
|
|||||||||||||||||||||||||
|
|
|
|
|
|
|
|
: : : |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
a0 n 1 + a1 n + + an 1 2n 2 + 2n 1 = 0; (j = n 1): |
||||||||||||||||||||||||||||||
|
Пример. [a; b] = [ 1; 1]; p(x) = 1; 0 = 2; 1 = 0; 2 = |
2 |
; n = 2 |
|
|
|
|
|
|
||||||||||||||||||||||
|
|
|
|
|
|
|
|
|
|||||||||||||||||||||||
|
|
3 |
|
|
|
|
|
|
|||||||||||||||||||||||
|
Тогда |
|
|
|
|
|
|
2 |
|
|
|
|
|
1 |
|
|
|
|
|
|
|
|
|
|
|
|
|
||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||
|
|
|
|
|
2a0 + |
|
|
= 0; a0 = |
|
; |
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||
|
|
|
|
|
|
3 |
3 |
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||||
|
|
|
|
|
2 |
a1 = 0; a1 |
= 0; |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||
|
|
|
|
|
|
3 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||
|
|
|
|
|
|
|
|
1 |
|
+ x2: |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||
|
|
|
|
|
|
2(x) = |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||
|
|
|
|
|
|
3 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||
|
Òàê êàê n(x) определяется однозначно, то эта система имеет единственное реше- |
||||||||||||||||||||||||||||||
íèå. |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
2. Находим n нулей построенной функции |
|
n(x) и эти нули возьмем в качестве |
||||||||||||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1 |
|
|
|
1 |
|
||
узлов интерполяционной квадратурной формулы. (x1 = p |
|
; |
x2 = p |
|
). Òàê êàê |
||||||||||||||||||||||||||
3 |
3 |
||||||||||||||||||||||||||||||
квадратурная формула точна для любого полинома степени |
(n 1), òî äëÿ Ak ïî- |
||||||||||||||||||||||||||||||
10
лучаем систему уравнений
A1 + A2 + + An = 0
A1x1 + A2x2 + + Anxn = 1
: : :
A1xn1 1 + A2xn2 1 + + Anxnn 1 = n
|
|
A1 + A2 = 2 |
|
|
|
||||||
|
|
1 |
|
1 |
|
|
|
|
|
||
|
|
A1( p |
|
) + A2(p |
|
) = 0 |
|
|
|
|
|
|
|
3 |
3 |
|
|
|
|||||
|
|
A1 = 1; A2 = 1: |
|
|
|
||||||
Квадратурная формула |
|
|
|
|
|
|
|
||||
|
b |
n |
1 |
|
|
|
|
|
|
||
Z |
|
Z f(x) dx f( p3) + f(p3); |
|||||||||
|
f(x) dx k=1 Akf(xk); |
||||||||||
|
|
X |
1 |
1 |
|
||||||
a |
|
1 |
|
|
|
|
|||||
построенная таким образом, точна для полиномов степени (2n 1). Действительно, любой полином Q2n 1(x) степени (2n 1) можно представить в виде
|
|
|
Q2n 1(x) = |
n(x)[полином степени (n 1)] + Pn 1(x) |
|||
è |
|
|
|
|
|
|
|
|
Z |
b |
|
n |
b |
n |
|
|
|
|
|
||||
|
p(x)Q2n 1 |
(x) dx = k=1 |
AkQ2n 1(xk) = Z p(x)Pn 1(x) dx k=1 AkPn 1(xk) = 0: |
||||
|
a |
|
|
X |
a |
X |
|
|
Среди полиномов степени 2n укажем полином Q2n(x) |
= 2(x), для которого |
|||||
b |
|
|
|
|
|
n |
|
p(x) n2(x) dx > 0, íî |
kn=1 Ak n2(xk) = 0 и стало быть (2n 1) наивысший по- |
||||||
a |
|||||||
R |
|
|
P |
|
|
||
рядок алгебраической точности квадратурной формулы типа Гаусса.
Отметим свойство коэффициентов Ak, важное при оценке неустранимой погрешности квадратурной формулы:
Ak > 0:
11
|
|
n(x) |
2 |
является полиномом степени (2n 2) и |
|
Действительно, функция f(x) = ( |
|
) |
|
||
x xk |
|
||||
равна нулю во всех узлах xi кроме узла xk è f(xk) > 0. |
|||||
b |
n |
|
|
||
Ra |
iP |
|
|
||
Тогда p(x)f(x) dx = |
Aif(xk) и значение Ak > 0. |
||||
|
=1 |
|
|
|
|
Методическая погрешность квадратурной формулы. (f 2 C2n[a; b])
Для функции f 2 C2n[a; b] построим интерполяционную таблицу
f(xk) = fk; f0(xk) = fk0 ; k = 1; 2; : : : ; n
и построим интерполяционный полином Эрмита H2n 1(x), удовлетворяющий этой таблице. Для этого полинома верна оценка
|
|
|
|
|
|
|
|
jf(x) H2n 1(x)j |
1 |
|
M2n n2 (x); n ! n: |
|
|
|
|
|
|
||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||||||
|
|
|
|
|
|
|
|
(2n)! |
|
|
|
|
|
|
|||||||||||||
|
|
|
|
b |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
b |
|
|
|
||
|
Òàê êàê |
|
p(x)H2n |
1(x) dx = |
n |
|
AkH2n |
|
1(xk) = |
n=1 Akf(xk), òî |
j Ra |
p(x)f(x) dx |
|
||||||||||||||
|
|
|
|
Ra |
|
|
|
b |
Pk=1 |
|
|
|
|
|
|
|
Pk b |
|
|
|
|
|
|||||
|
|
n |
A |
f(x ) |
|
= |
p(x)[f(x) |
|
|
H |
|
|
(x)] dx |
j Ra |
p(x) f(x) |
|
H |
|
(x) |
dx |
|
||||||
|
|
k=1 |
k |
|
b k |
j |
|
j Ra |
|
|
|
2n 1 |
|
|
|
j |
|
|
2n 1 |
j |
|
|
|||||
|
P1 |
|
|
Ra p(x) |
n2 (x) dx; ! |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||
|
M2n |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||
(2n)! |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||
Замечание. Ортогональные полиномы нами были введены в главе 6 (Аппроксимация функций, 2). Там же приведены примеры ортогональных полиномов (полиномы Лежандра, Чебышева, Лаггера, Эрмита). Ортогональные полиномы образуют базисы в сепарабельных гильбертовых пространствах функций. Кроме того они являются ненулевыми решениями однородных линейных дифференциальных уравнений второго порядка.
Сходимость квадратурного процесса Гаусса.
Для заданной функции f 2 C[a; b] рассмотрим последовательность квадратурных
формул Гаусса и величин методических погрешностей Rn(f)
Rn(f) = Z |
b |
n |
|
|
|
||
|
p(x)f(x) dx k=1 |
Ak(n)f(xk); |
|
a |
|
X |
|
ãäå x1; x2; : : : ; xn нули ортогональных полиномов n.
Теорема. Для любой заданной непрерывной функции f(x) значения jRn(f)j ! 0
12
ïðè n ! 1:
Доказательство. Для любого " > 0 для достаточно больших n можно найти поли-
íîì P2n 1(x) степени (2n 1) такой, что
" 1 jf(x) P2n 1(x)j < 2 Rb
p(x) dx
a
для всех x 2 [a; b]: Зафиксировав n, построим квадратурную формулу типа Гаусса. Для нее
|
b |
|
|
|
|
b |
|
|
|
|
|
b |
|
|
|
|
n |
||
jRn(f)j j Z p(x)f(x) dx Z |
p(x)P2n 1(x) dxj+j Z |
|
|
|
|
||||||||||||||
p(x)P2n 1(x) dx k=1 Ak(n)P2n 1(xk)j+ |
|||||||||||||||||||
a |
|
|
|
|
|
a |
|
|
|
|
|
a |
|
|
|
|
X |
||
n |
|
|
|
|
|
|
|
|
n |
|
|
|
|
|
n |
|
|
|
|
X |
(n) |
|
|
|
|
|
|
Xk |
(n) |
|
|
" |
|
X |
(n) |
|
|
|
|
+j |
|
Ak |
P2n 1(xk) |
|
|
Ak |
f(xk)j |
|
+ |
|
Ak |
jP2n 1 |
(xk) |
f(xk)j |
|||||
k=1 |
|
|
|
|
|
|
|
=1 |
|
|
|
|
|
k=1 |
|
|
|
|
|
|
|
|
" |
|
" |
|
n |
(n) |
1 |
|
|
|
(n) |
|
|
||||
|
|
|
|
|
|
|
|
|
Xk |
|
|
|
|
|
|
|
|
|
|
|
|
|
2 + 2 |
Ak |
|
p(x) dx |
(òàê êàê Ak > 0): |
|
|||||||||||
|
|
|
|
|
b |
|
|||||||||||||
|
|
|
|
|
|
|
|
|
=1 |
|
Ra |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
b |
p(x) 1 dx = Pkn=1 Ak(n), òî jRn(f)j < " при достаточно больших n. |
||||||||||||||||||
Òàê êàê Ra |
|||||||||||||||||||
Таким образом, сходимость квадратурного процесса основана на положительности коэффициентов Ak. Положительность этих коэффициентов влечет и устойчивость квадратурных формул типа Гаусса относительно погрешностей значений функции f(x) в узлах xk: åñëè j f(xk)j < k < , то неустранимая погрешность квадратурной формулы оценивается величиной
n |
n |
|
b |
|
= Z p(x) dx |
||
k=1 |
Ak(n)j kj k=1 |
Ak(n) |
|
X |
X |
a |
|
и не возрастает с ростом n.
13
7.3Апостериорная оценка погрешности квадратурных формул (правило Рунге)
Можно считать (как мы видели во Введении), что погрешность Rh(f) квадратурной формулы с шагом h имеет вид (при малых h)
Rh(f) = Cmhm + o(hm);
где число m определяет число непрерывных производных функции f, а величина Cm слабо зависит от h.
Обозначим Jh значение интеграла
b
Z
J = p(x)f(x) dx;
a
вычисленного приближенно по квадратурной формуле с шагом h,
J = Jh + Cmhm + o(hm):
Тогда J = Jh + Cm( |
h |
)m + o(hm): |
|
2 |
|||
2 |
|
Зная вычисленные значения Jh è Jh можно найти приближенное значение Cm,
|
|
1 |
2 |
|
|
|
|
|
|
|
|
|
|
|
пренебрегая величинами |
|
o(hm): |
|
|
|
|
|
|
|
|
|
|
||
|
m |
|
|
|
|
|
|
|
|
|
|
|||
|
|
h |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
J |
2 |
|
Jh |
|
||||||
|
|
|
|
|
h |
|
||||||||
|
|
|
Cm = |
|
|
|
|
|
|
|
: |
|
||
|
|
|
hm(1 |
1 |
) |
|
||||||||
|
|
|
|
2m |
|
|
|
|||||||
Тогда |
|
|
|
|
|
J 2 Jh |
m |
|||||||
|
|
|
|
|
|
|||||||||
|
|
|
|
|
|
|
h |
|
||||||
|
Rh(f) = J Jh = |
|
|
|
|
|
|
+ o(h ) |
||||||
|
|
|
1 |
1 |
|
|||||||||
|
|
2m |
|
|||||||||||
è |
|
h |
h |
|
|
J 2 Jh |
m |
|||||||
|
|
|
|
|||||||||||
|
|
|
|
|
|
|
h |
|
||||||
|
R 2 (f) = J J 2 = |
|
+ o(h ) |
|||||||||||
|
2m 1 |
|||||||||||||
и получаем апостериорные оценки значений Rh(f) è Rh (f) по правилу Рунге.
2
Деля шаг h на L частей, получаем аналогичным способом априорные оценки погрешностей Rh(f) по правилу Рунге.
Правило Рунге используется для приближенного вычисления интегралов с авто-
14
матическим выбором шага. Для каждого частичного интервала [ai; bi] погрешности оцениваются по правилу Рунге при дроблении исходного шага на этом интервале.
Существуют и другие варианты апостериорных оценок погрешностей квадратурных формул, основанные на других представлениях Rh(f) (методы Ричардсона, см. ¾Методические указания по вычислительному практикуму¿, И. В. Олемской).
15
Глава 8 Методы решения задачи Коши для
систем обыкновенных дифференциальных уравнений
ВВЕДЕНИЕ.
Рассматривается задача Коши для системы обыковенных дифференциальных урав-
нений для функций ui(t) du (t)
8 |
i |
= fi(t;0u1(t); u2(t); : : : ; un(t)); t 2 [0; T ]; |
||
dt |
||||
<ui(t = 0) = ui : |
||||
: |
|
|
|
|
Обозначая u(t) = fu1(t); : : : ; un(t)g, запишем систему в виде |
||||
|
du |
|||
|
|
|
= f(t; u(t)); t 2 [0; T ]; u(0) = u0: |
|
|
|
dt |
||
Большинство численных методов решения этой задачи можно получить, исходя из идей решения задачи Коши для скалярной функции u(t).
Рассмотрим эквивалентное нелинейное интегральное уравнение
u(t) = u0 + Z0 |
t |
f( ; u( )) d : |
Из курса дифференциальных уравнений известно, что при условиях, что функция
16
f(t; u) непрерывна как функция двух переменных и удовлетворяет условию Липшица jf(t; u00) f(t; u0)j < Lju00 u0j; решение задачи Коши существует, единственно и последовательность fuk(t)g
t
Z
uk(t) = u0 + f( ; uk 1( )) d ; u0(t) = u0
0
равномерно сходится к точному решению u (t) (теорема Пикара).
Можно строить численные методы для получения приближенных значений uk(t) в узлах ti 2 [0; T ] на основе квадратурных формул.
Можно также строить методы, исходя из формул численного дифференцирования, рассматривая узлы ftig на отрезке [0; T ].
На простом примере введем два важных понятия.
Пример. Рассмотрим разностную схему (явная разностная схема Эйлера):
8ui+1 ui
< = f(ti; ui); ti = i; i = 0; 1; : : : ; N 1; N = T;
:u0 = u0:
Значения ui+1 получаем последовательно |
ui+1 = ui + f(ti; ui). Обозначим ui |
||||||||||||||||||
значения решения дифференциальной задачи в узлах ti |
и рассмотрим разность |
||||||||||||||||||
ui ui = zi. |
|
|
|
|
du (ti + ) |
|
|
|
|
|
|
|
|
|
|
|
|||
Òàê êàê u |
|
= u + |
|
|
= u |
+ f(t |
|
+ )); |
0 |
6 |
|
6 |
1, а функция |
||||||
|
|
dt |
|
i |
|||||||||||||||
|
i+1 |
i |
|
|
|
|
|
i |
|
|
|
|
|
||||||
f(ti; ui ) непрерывна по t, то |
|
|
|
|
|
|
|
|
|
|
|
||||||||
|
|
|
|
|
|
|
|
ui+1 ui |
= f(ti; u ) + O( ): |
|
|
|
|
|
|||||
|
|
|
|
|
|
|
|
|
|
|
|
i |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Обозначим |
|
= f(t |
; u ) |
|
ui+1 ui |
. |
|
|
|
|
|
|
|
|
|
||||
|
i |
|
i |
i |
|
|
|
|
|
|
|
|
|
|
|
|
|
||
Эта величина определяет погрешность разностной схемы в узлах ti для точного решения u (t). В нашем примере величина i = O( ).
Говорят, что порядок аппроксимации разностной схемы Эйлера равен 1 èëè разностная схема имеет первый порядок аппроксимации (первый порядок точности).
17
Рассмотрим сеточную функцию zi = ui ui. ßñíî, ÷òî
zi+1 zi |
= |
ui+1 ui |
|
ui+1 ui |
= f(t |
; u ) |
|
i |
f(t |
; u |
) |
|
|
|
|
|
|
i |
i |
i |
i |
|
|||
è zi+1 zi = [f(ti; ui ) f(ti; ui)] |
i. Согласно предположению относительно |
|||||||||||
свойств функции f(t; u) |
|
|
|
|
|
|
|
|
|
|
||
jf(ti; ui ) f(ti; ui)j 6 Ljui uij = Ljzij:
Тогда jzi+1j 6 jzij + Ljzij + j ij èëè jzi+1j 6 (1 + L)jzij + j ij. Далее
jzi+1 6 (1+ L)2jzi 1j+ [(1+ L)j i 1j+j ij] 6 (1+ L)2jzi 1j+2(1+ L) maxfj ij; j i 1jg
(грубые оценки).
Продолжая этот процесс, получаем
|
|
z |
i+1j 6 (1 + |
L |
i+1 |
z |
0j + (1 + |
i |
|
(1 + |
L i |
|
max |
j |
kj |
6 |
T (1 + L)i |
max |
; |
|||||
|
|
j |
) |
|
|
j |
|
) |
|
) |
|
k |
|
k |
j ij |
|||||||||
è òàê êàê z0 = 0, òî |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||
|
|
|
|
|
|
|
|
|
z |
i+1j 6 |
T |
(1 + |
L |
i |
|
max |
j |
kj |
: |
|
|
|
||
|
|
|
|
|
|
|
|
|
j |
|
|
) |
|
k |
|
|
|
|
||||||
|
|
Оценим порядок величины (1 + L)i = y при ! 0: ln y = i ln(1 + L) и для малых |
||||||||||||||||||||||
|
: |
|
~ |
, тогда |
y e |
T L. |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||
|
ln(y)i L |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||
|
|
Таким образом jui uij 6 TeLT maxk j kj |
|
= TeLT O( ) è ïðè |
! 0 значения |
|||||||||||||||||||
maxi jui uij ! 0: |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||
Устойчивость разностной схемы относительно вариации правых частей f(t; u). Наряду с разностной схемой для fuig рассмотрим разностную схему
8vi+1 vi
< = f(ti; vi) + "i;
:v0 = u0:
Относительно значений "i будем предполагать, что j"ij < ". Обозначим разность решений vi ui = wi. Для сеточной функции получим
wi+1 wi = f(ti; vi) f(ti; ui) "i:
18
Òàê êàê jf(ti; vi) f(ti; ui)j 6 Ljvi uij = Ljwij, òî jwi+1j 6 (1 + L)jwij + j"ij.
Повторяя процедуры предыдущего пункта, получим maxi jvi uij 6 T eT L", и разностная явная схема Эйлера устойчива относительно вариаций правых частей:
max jvi uij = O("):
i
Теорема Филлипова. Из аппроксимации разностной схемы и ее устойчивости следует сходимость разностной схемы к решению fui g.
Действительно, jui vij 6 jui uij+ jui vij è maxi jvi ui j = T eT L(O( ) + O(")) = = O( ) + O(").
8.1Одношаговые методы (методы Рунге-Кутты)
Усложнив разностную схему, можно повысить порядок аппроксимации. Исходя из явной разностной схемы Эйлера, значение un+1 будем находить в два этапа:
I этап: схема Эйлера применяется для шага :
u~n = un + f(tn; un); 0 < < 1:
II этап: схема Эйлера применяется для шага , но значение функции f вычисляется
ïðè t = tn è ïðè t = tn + и берется их линейная комбинация с весам 1 и :
un+1 = un + [(1 )f(tn; un) + f(tn + ; u~n)]
В результате получим разностную схему:
un+1 un = (1 )f(tn; un) + f(tn + ; un + f(tn; un)):
Порядок аппроксимации этой разностной схемы определим, оценивая величину
n: |
|
|
|
|
|
|
|
|
un+1 un |
|
n = (1 |
|
f t |
; u |
f t |
; u |
+ |
f t |
; u |
: |
|
) |
( n |
n) + |
( n + |
n |
( n |
n)) |
|
|
Предполагая, что творая проиводная решения u (t) непрерывна, получаем по фор-
19