1курс,2семестр лабы для зачета / Метода к лабам
.pdfy y2 х, y(1) 0
на отрезке x [1,3] с шагом по сетке h 0.2 .
Из постановки задачи видно, что f (x, y) y y2 х , а из формулы (12.8)
следует, что
yk 1 yk h f (x, y) .
Зная, что y(1) 0, x0 1, получаем первое значение: y1 y0 h f (x0 , y0 ) 0 0.2 (( y0 )2 x0 ) 0.2.
Далее продолжаем процесс рассчета по формулам
x1 1.2, y1 0.2 y2 y1 h f (x1, y1) 0.2 0.2 (( y1)2 x1) 0.432; x2 1.4, y2 0.432 y3 y2 h f (x2 , y2 )
0.432 0.2 (( y2 )2 x2 ) 0.6746752;
...
x10 2.8, y10 1.56185 y11 y10 h f (x10 , y10 )
1.56185 0.2 (( y10 )2 x10 ) 1.634; |
|
|||
x11 3, y11 1.634. |
|
|
|
|
М2. Неявная схема 1-го порядка |
|
|
||
Используя в (12.5) формулу правых прямоугольников, получим |
|
|||
|
yk 1 yk |
f (x |
, yk 1). |
(12.9) |
|
|
|||
|
h |
k 1 |
|
|
|
|
|
|
Эта схема неразрешима явно относительно yk 1, поэтому для получения требуется использовать итерационную процедуру решения уравнения
(12.9) (см. метод простой итерации в подразд. 9.2):
yk 1,s yk f (xk 1, yk 1,s 1),
где s 1, 2, 3, ... – номер итерации.
За начальное приближение берется значение yk 1, 0 yk c предыдущего шага. Обычно, если значение h выбрано удачно, достаточно сделать 2–3 итера-
ции для достижения заданной погрешности |
|
|
. Эффектив- |
|
|
yk 1, s yk 1, s 1 |
|
ность неявной схемы заключается в том, что у нее константа устойчивости С0 значительно меньше, чем у явной схемы.
M3. Неявная схема второго порядка
Используя в (12.5) формулу трапеций, получим
yk 1 yk |
f (x, yk ) f (x |
, yk 1) |
|
|
|
|
|
k 1 |
|
. |
(12.10) |
|
|
|
|||
h |
2 |
|
|
|
101
Т. к. формула трапеций имеет третий порядок точности на интервалеxx , xk 1 , то погрешность аппроксимации (h) – второй.
Схема (12.10) не разрешена относительно yk 1, поэтому требуется итерационная процедура (см. схему М2):
yk 1, s yk |
h |
( f (x , |
yk ) f (x |
|
, yk 1, s 1), |
s 1, 2, ..., |
yk 1, 0 |
yk . |
|
k 1 |
|||||||
2 |
k |
|
|
|
|
|
||
|
|
|
|
|
|
|
Алгоритм представлен на рис. 12.2.
Рис. 12.2
Рассмотрим пример решения задачи Коши неявным методом Эйлера второго порядка (с пересчетом).
Задача Коши имеет вид:
102
|
y y2 х, |
y(1) 0 |
||
на отрезке x [1, 3] с шагом по сетке h 0.2 . |
||||
Из постановки задачи |
видно, что |
f (x, y) y y2 х , а из формулы |
||
(12.10) следует, что для начала необходимо рассчитать прогноз |
||||
|
yi 1 yi h f (xi , yi ) |
|||
а затем произвести коррекцию: |
|
|
||
y |
y h f (xi , yi ) f (xi 1, yi 1) . |
|||
i 1 |
i |
|
2 |
|
|
|
|
|
Рассчитаем значение по схеме для первой точки отрезка по х:
x0 1, y0 0; y1 |
y0 |
h f (x0 , y0 ) 0 0.2 (0 0 1) 0.2; |
|
|
|||||||||||||
|
|
|
|
|
|
f (x , y ) f (x , y ) |
|
|
|
(0 0 1) (( 0.2)2 |
1.2) |
|
|||||
y |
|
y |
h |
|
|
|
0 |
0 |
1 1 |
|
0 0.2 |
|
|
|
0.216; |
||
|
|
|
|
|
|
|
|
|
|
||||||||
1 |
0 |
|
|
|
|
|
|
|
2 |
|
|
|
2 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
x1 1.2, y1 0.216; |
|
|
|
|
|
|
|
|
|||||||||
y |
2 |
y |
h f (x , y ) 0.216 0.2 (( 0.216)2 1.2) 0.44667; |
|
|||||||||||||
|
1 |
|
|
1 |
1 |
|
|
|
|
|
|
|
|
||||
y |
2 |
y h |
f (x1, y1) f (x2 , y2 ) |
|
|
|
|
|
|
||||||||
|
|
|
|
|
|
|
|||||||||||
|
1 |
|
|
|
|
|
|
|
2 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
0.216 0.2 |
(( 0.216)2 1.2) (( 0.44667)2 1.4) |
0.452; |
|
||||||||||||
|
|
|
|
2 |
|
|
|
||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
Продолжая, мы приходим к значению в последней точке |
|
|
|||||||||||||||
|
x10 2.8, y10 1,54075; |
|
|
|
|
|
|
|
|||||||||
|
y y |
h f (x |
, y ) 1,54075 0.2 (( 1,54075)2 2.8) 1,62597; |
||||||||||||||
|
|
11 |
10 |
|
|
|
|
|
10 |
10 |
|
|
|
|
|
|
|
|
y y |
h |
f (x10 , y10 ) f (x11, y11) |
|
|
|
|
|
|||||||||
|
|
|
|
|
|
||||||||||||
|
|
11 |
10 |
|
|
|
|
|
|
2 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1,54075 0.2 (( 1,54075)2 2.8) (( 1,62597)2 3) 1.61898. 2
Легко заметить, что данные результаты схожи с результатами, полученными с помощью явного метода Эйлера по схеме М1.
M4. Схема предиктор – корректор (Рунге – Кутта) второго порядка
Используя в (12.5) формулу средних, получим |
|
|
|||
|
yk 1 yk |
f (x |
, yk 1/ 2 ) . |
(12.11) |
|
|
|
||||
|
h |
k 1/ 2 |
|
|
|
|
|
|
|
|
|
Уравнение разрешено явно относительно yk 1, однако в правой части |
|||||
присутствует неизвестное значение yk 1/ 2 |
в середине отрезка [x , |
x |
]. Для |
||
|
|
|
k |
k 1 |
|
решения этого уравнения существует следующий способ. Вначале по явной схеме (12.8) рассчитывают yk 1/ 2 (предиктор):
103
yk 1/ 2 yk h2 f (xk , yk ) .
После этого рассчитывают yk 1 по (12.11) (корректор). В результате схе-
ма оказывается явной и имеет второй порядок. Заметим, что схема получается из схемы М3, если в ней выполнять только две итерации (s = 1, 2).
Алгоритм представлен на рис. 12.3.
Начало |
a, b, nx, np, ny, |
y, FPR, OUT |
h=(b – a)/nx |
OUT x, y, ny |
n=1, nx |
FPR x, y, f pr |
n=1, ny |
y pri yi h 2 f pri |
x+=h/2 |
FPR x, y pr , f , ny |
i=1, ny |
yi=yi+h fi |
n mod np=0 |
OUT x, y, np |
x+=h/2 |
Конец |
Рис. 12.3 |
М5. Схема Рунге – Кутта четвертого порядка |
|
|
|
||||||
Используя в (12.5) формулу Симпсона, получим |
|
|
|
||||||
|
yk 1 yk |
|
1 |
[ f (x , |
yk ) 4 f (x |
, yk 1/ 2 ) f (x |
, |
yk 1)] . |
(12.12) |
|
|
|
|||||||
|
h |
|
6 |
k |
k 1/ 2 |
k 1 |
|
|
|
|
|
|
|
|
|
|
|
Т. к. формула Симпсона имеет пятый порядок точности, то погрешность аппроксимации данной схемы имеет четвертый порядок.
104
Можно по-разному реализовать расчет неявного по yk 1 уравнения
(12.12), однако наибольшее распространение получил следующий способ. Делают предиктор вида
k0 f (xi , y); |
|
|
|
|||
k |
|
yk 1/ 2,1 |
yk h / 2 f (x , y); |
|
||
1 |
|
|
|
i |
|
|
k |
2 |
yk 1/ 2,2 |
yk h / 2 |
f (x |
, yk 1/ 2,1); |
|
|
|
|
|
i 1/ 2 |
|
|
k |
3 |
yk 1,1 yk h f (x |
|
, yk 1/ 2,2 ), |
||
|
|
i 1/ 2 |
|
|
затем корректор по формуле
yk 1 yk |
h |
[ f (x , |
yk ) 2 f (x |
|
yk 1/ 2,1) |
|||
|
k 1/ 2, |
|||||||
|
6 |
k |
|
|
|
|
||
|
|
|
|
|
|
|
||
2 f (x |
k 1/ 2, |
yk 1/ 2, 2 ) f (x |
k 1 |
, yk 1, 1)]. |
||||
|
|
|
|
|
|
Алгоритм метода представлен на рис. 12.4.
|
Начало |
|
|
A |
|
B |
|
|
C |
|
|
|
|
|
|
|
|||
a, b, nx, np, ny, |
|
|
FPR x, ypr , k2 |
|
|
||||
|
|
|
|
|
|
|
|
||
y, FPR, OUT |
|
|
|
|
|
|
|
|
|
|
h=(b – a)/nx |
|
|
|
|
n=1, ny |
|
|
|
|
|
|
|
|
|
|
|
|
|
OUT x, y, ny |
|
|
y pr |
|
yi h 2 |
k2 |
|
|
|
|
|
|
|
|
i |
|
|||
|
|
|
|
i |
|
|
|
||
|
n=1, nx |
|
|
|
|
|
|
|
|
FPR x, y,k0 |
|
|
|
|
x+=h/2 |
|
|
|
|
|
|
|
|
|
|
|
|
||
|
n=1, ny |
|
|
FPR x, y pr , k2 |
|
||||
|
|
|
|
|
|
n=1, ny |
|
|
|
y pr |
yi h 2 k0 |
i |
|
|
|
|
|
|
|
i |
|
yi yi h 6 |
(k0 i 2 k1i |
2 k2 i |
k3 i ) |
||||
|
|
|
|
||||||
|
x+=h/2 |
|
|
|
|
|
|
|
|
FPR x, ypr ,k1 |
|
|
OUT x, y, np |
|
|
||||
|
|
|
|
|
|
|
|
||
|
n=1, ny |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Конец |
|
|
|
y pr |
yi h 2 k1 |
i |
|
|
|
|
|
|
|
i |
|
|
|
|
|
|
|
|
x+=h/2
A B C
Рис. 12.4
12.4. Многошаговые схемы Адамса
При построении всех предыдущих схем для вычисления интеграла в правой части (12.4) использовались лишь точки в диапазоне одного шага [xk , xk 1].
Поэтому при реализации таких схем для вычисления следующего значения yk 1
105
требуется знать только одно предыдущее значение yk , т. е. рекуррентная последовательность получается первого порядка. Такие схемы называют одношаговыми. Мы, однако, видели, что для повышения точности при переходе от xk к xk+1 приходилось использовать и значения функции F внутри интервала [xk 1/ 2, yk 1/ 2 ]. Схемы, в которых это используется (M4, M5, ...), называют схе-
мами с дробными шагами. В этих схемах повышение точности достигается за счет дополнительных затрат на вычисление функции F(x) в промежуточных точках интервала [xk , xk 1].
Идея методов Адамса заключается в том, чтобы для повышения точности использовать уже вычисленные на предыдущих шагах значения yk , yk 1, yk 2, ... .
Заменим в (12.4) F(x) интерполяционным многочленом Ньютона вида
F (x) F (xk ) (x xk ) F (xk ) hF (xk 1)
(x xk )(x xk 1) F (xk ) 2F (xk 1) F (xk 2 ) ... .
2h2
После интегрирования на интервале [xk , xk 1] получим явную экстрапо-
ляционную схему Адамса (экстраполяцией называется получение значений интерполяционного многочлена в точках x, выходящих за крайние узлы сетки). В нашем случае интегрирование производится на интервале [xk , xk 1], а поли-
ном строится по узлам xk , xk 1, xk 2 .
Порядок аппроксимации схемы в этом случае определяется количеством использованных при построении полинома узлов (например если используются xk , xk 1, то схема второго порядка).
Если в (12.4) F(x) заменить многочленом Ньютона вида
F (x) F (xk 1) (x xk 1) F (xk 1)h F (xk )
(x xk 1)(x xk ) F (xk 1) 2F (xk ) F (xk 1) ... , 2h2
то после интегрирования получим неявную интерполяционную схему Адамса.
Заметим, что неявная интерполяционная схема второго порядка совпадает со схемой М3.
M6. Явная экстраполяционная схема Адамса второго порядка
yk 1 yk |
1.5 f (x , |
yk ) 0.5 f (x |
, yk 1) . |
(12.13) |
|
||||
h |
k |
k 1 |
|
|
|
|
|
|
Схема двухшаговая, поэтому для начала расчетов необходимо найти y1 по методу М4, после чего y2 , y3, ... вычислять по (12.13).
106
М7. Явная экстраполяционная схема Адамса третьего порядка
|
yk 1 yk |
|
23 |
|
f (x , |
yk ) |
16 |
|
|
f (x |
, yk 1) |
5 |
f (x |
|
, yk 2 ) . (12.14) |
||||
|
|
|
|
|
|
|
|
||||||||||||
|
h |
12 |
|
k |
12 |
|
|
|
k 1 |
|
12 |
k 2 |
|
|
|||||
|
|
|
|
|
|
|
|
|
|
|
|
||||||||
Схема трехшаговая, поэтому для |
|
начала расчетов необходимо |
найти |
||||||||||||||||
y1, y2 по методу М5, после чего y3, y4, ... вычислять по (12.14). |
|
||||||||||||||||||
M8. Неявная схема Адамса третьего порядка |
|
|
|
|
|
|
|||||||||||||
|
yk 1 yk |
|
5 |
|
f (x |
, yk 1) |
|
8 |
f (x , yk ) |
1 |
f (x |
, yk 1) . (12.15) |
|||||||
|
|
|
|
|
|
||||||||||||||
|
h |
12 |
|
|
k 1 |
|
|
|
12 |
|
k |
12 |
k 1 |
|
|
|
|||
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||
Т. к.схема двухшаговая, то для начала расчетов необходимо найти |
y1 по |
||||||||||||||||||
методу М5, после чего y2 , |
y3, ... вычислять по (12.15). |
|
|
|
|
|
|
Для нахождения yk 1 требуется использовать метод простой итерации:
yk 1, s
Значение
yk h |
5 |
f (x |
, yk 1, s 1) |
|
8 |
f (x , yk ) |
1 |
|
f (x |
|
|
, yk 1) |
. |
|||||||||||||
|
|
|
|
|
|
1 |
||||||||||||||||||||
|
|
|
|
k 1 |
|
|
12 |
k |
12 |
|
|
k |
|
|
||||||||||||
|
|
|
12 |
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||
yk 1, 0 |
следует рассчитать по формуле (12.13): |
|
|
|
|
|
|
|
||||||||||||||||||
y |
k 1, 0 |
y |
k |
h |
|
f (x , |
y |
k |
) 0.5 |
f (x |
, |
|
y |
k 1 |
) |
|
. |
|
||||||||
|
|
|
|
1.5 |
|
|
|
|
|
|||||||||||||||||
|
|
|
|
|
|
|
|
|
|
k |
|
|
|
|
|
k 1 |
|
|
|
|
|
|
|
|
Чаще всего бывает достаточно одной итерации. Если при этом разность yk 1,0 yk 1,1 оказывается большой, то следует уменьшить h. Схема алгоритма представлена на рис. 12.5.
107
Рис. 12.5. 1-й фрагмент (окончание см. на с. 109)
108
Рис. 12.5 окончание (начало см. на с. 108)
109
12.5. Особенности программирования алгоритмов
Для решения задачи Коши составляется стандартная подпрограмма, которая помещается в библиотеку стандартных программ. Передача всех необходимых данных из основной программы пользователя в подпрограмму организуется через список формальных параметров, который включает:
a, b – начало и конец области интегрирования; nx – количество шагов сетки;
np – параметр, указывающий, через какое количество шагов организовывать вывод данных;
ny – количество уравнений;
y – массив, в который сначала при обращении к подпрограмме помещается вектор начальных значений y(0) u0 , а по окончании работы в нем помещается решение на конце интервала y(b).
Для итерационных методов дополнительно вводятся: nit – максимально допустимое количество итераций;
– точность (погрешность, до которой выполняются итерации).
В подпрограмме FPR(x, y, F, ny) для заданных x и y вычисляется вектор f f (x, y), в подпрограмме OUT (x, y, ny) осуществляется вывод данных.
12.6. Индивидуальные задания
Составить отдельную подпрограмму, оформленную в виде модуля, для решения задачи Коши в соответствии со схемой согласно варианту (табл. 12.1).
С помощью этой подпрограммы решить задачу для системы двух уравнений в соответствии с вариантом (см. табл. 12.1):
dudx1 f1(x, u1, u2 ); dudx2 f2 (x, u1, u2 ); a x b;
u1(a) u10; u2 (a) u20.
Точное решение этой задачи при u10 2a, u20 ea одинаково для всех вариантов и имеет вид u1 2x, u2 ex .
Функция FPR для варианта 15 имеет вид
void FPR(double x, double *y,double *f)
{
f[0]=y[0]*exp(x)/(x*y[1]); f[1]=2*x/y[0]+y[1]-1;}
Вариант процедуры OUT, обеспечивающей вывод таблицы значений приближенного и точного решения и погрешностей ( d1, d2 ), имеет вид
110