ВУЗ: Не указан

Категория: Не указан

Дисциплина: Не указана

Добавлен: 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