ВУЗ: Не указан
Категория: Не указан
Дисциплина: Не указана
Добавлен: 10.04.2024
Просмотров: 82
Скачиваний: 0
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