Математична основа способу Мілна

1 2 3 4 5 6

Вихідні дані для реалізації системи звичайних диференційних рівнянь

Система диференційних рівнянь:

Початкові умови: А=0, В=1

t(0)=0, x(0)=0, y(0)=0

Задана точність: Е=

Обираємо с=6

Математична основа засобу Рунге–Кутта

Засіб Рунге -Кутта можливо получити, якщо разкласти у ряд Тейлора значення у(х)

y(x₀+h)=y(x₀)+h(x₀)h³ +hⁿyⁿ(x₀)

x_i=x(0)+Ih

y_i+1=y_i+∙(K_1i+2K_2i+2K_3i+2K_4i)

K_1i=h∙f(x_i,y_i)

K_2i=h∙f(x_i+∙y_i+)

K_3i=h∙f(x_i+

K_4i=h∙f(x_i+h∙y_i+K₃)

Блок - схема головного модуля по Рунге - Кутту:

Реалізація програми за засобом Рунге – Кутта.:

DECLARE SUB KUTT (T!, X!, Y!, A%, B%, C%, E!, H!, N%, T(), X(), Y(), K1X!, K1Y!, K2X!, K2Y!, K3X!, K3Y!, K4X!, K4Y!)

DECLARE SUB GRAF (T!, X!, Y!, A%, B%, C%, E!, H!, N%, T(), X(), Y(), K1X!, K1Y!, K2X!, K2Y!, K3X!, K3Y!, K4X!, K4Y!)

INPUT "C"; C%

E! = C% * 10 ^ (-4)

H! = E! ^ (1 / 4)

CONST A% = 0: CONST B% = 1

DIM SHARED T!(2000), X!(2000), Y!(2000), K1X!(2000), K1Y!(2000), K2X!(2000), K2Y!(2000), K3X!(2000), K3Y!(2000), K4X!(2000), K4Y!(2000)

T(0) = 0: X(0) = 0: Y(0) = 0

M1: CALL KUTT(T!, X!, Y!, A%, B%, C%, E!, H!, N%, T(), X(), Y(), K1X!, K1Y!, K2X!, K2Y!, K3X!, K3Y!, K4X!, K4Y!)

FOR I% = 0 TO N%

X1(I%) = X(I%)

Y1(I%) = Y(I%)

NEXT I%

H! = H! / 2

CALL KUTT(T!, X!, Y!, A%, B%, C%, E!, H!, N%, T(), X(), Y(), K1X!, K1Y!, K2X!, K2Y!, K3X!, K3Y!, K4X!, K4Y!)

FOR I% = 0 TO N%

IF ABS(X1(I%) - X(I%)) * (16 / 15) > E! THEN

GOTO M1

ELSE GOTO M2

END IF

IF ABS(Y1(I%) - Y(I%)) * (16 / 15) > E! THEN

GOTO M1

ELSE GOTO M2:

END IF

NEXT I%

M2: FOR I% = 1 TO N%

PRINT T(I%), X(I%), Y(I%)

NEXT I%

PRINT "H"; H!

INPUT K!

CALL GRAF(T!, X!, Y!, A%, B%, C%, E!, H!, N%, T(), X(), Y(), K1X!, K1Y!, K2X!, K2Y!, K3X!, K3Y!, K4X!, K4Y!)

END

SUB GRAF (T!, X!, Y!, A%, B%, C%, E!, H!, N%, T(), X(), X1(), K1X!, K1Y!, K2X!, K2Y!, K3X!, K3Y!, K4X!, K4Y!)

SCREEN 2

VIEW (170, 50)-(470, 150)

WINDOW (-1, 1.5)-(1, -1.5)

FOR I% = 0 TO N% - 1

PSET (T(I%), X(I%))

PSET (T(I%), Y(I%))

LINE (T(I%), X(I%))-(T(I% + 1), X(I% + 1))

LINE (T(I%), Y(I%))-(T(I% + 1), Y(I% + 1))

NEXT I%

LINE (-1, 0)-(1, 0)

LINE (0, -1.5)-(0, 1.5)

END SUB

SUB KUTT (T!, X!, Y!, A%, B%, C%, E!, H!, N%, T(), X(), Y(), K1X!, K1Y!, K2X!, K2Y!, K3X!, K3Y!, K4X!, K4Y!)

N% = (B% - A%) / H!

FOR I% = 0 TO N%

T(I%) = T(0) + I% * H!

K1X(I%) = H! * (-2 * X(I%) + 5 * Y(I%))

K1Y(I%) = H! * ((EXP(.5 * Y(I%) + T(I%)) - EXP(-.5 * Y(I%) + T(I%))) / 3 +.5 * Y(I%))

K2X(I%) = H! * (-2 * (X(I%) + K1X(I%) / 2) + 5 * (Y(I%) + K1Y(I%) / 2))

K2Y(I%) = H! * ((EXP(.5 * (Y(I%) + K1Y(I%) / 2) + (T(I%) + H! / 2) - EXP(-.5 * Y(I%) + K1Y(I%) / 2) - (T(I%) + H! / 2))) / 3 +.5 * (Y(I%) + K1Y(I%) / 2))

K3X(I%) = H! * (-2 * (X(I%) + K2X(I%) / 2) + 5 * (Y(I%) + K2Y(I%) / 2))

K3Y(I%) = H! * ((EXP(.5 * (Y(I%) + K2Y(I%) / 2) + (T(I%) + H! / 2) - EXP(-.5 * Y(I%) + K2Y(I%) / 2) - (T(I%) + H! / 2))) / 3 +.5 * (Y(I%) + K2Y(I%) / 2))

K4X(I%) = H! * (-2 * (X(I%) + K3X(I%)) + 5 * (Y(I%) + K3Y(I%)))

K4Y(I%) = H! * ((EXP(.5 * (Y(I%) + K3Y(I%)) + (T(I%) + H!) - EXP(-.5 * Y(I%) + K3Y(I%)) - (T(I%) + H!))) / 3 +.5 * (Y(I%) + K3Y(I%) / 2))

X(I% + 1) = X(I%) + 1 / 6 * (K1X(I%) + 2 * K2X(I%) + 2 * K3X(I%) + K4X(I%))

Y(I% + 1) = Y(I%) + 1 / 6 * (K1Y(I%) + 2 * K2Y(I%) + 2 * K3Y(I%) + K4Y(I%))

NEXT I%

END SUB

Результати реалізації системи диференційних рівнянь за засобом Рунге – Кутта.

T	X	Y
0	0	0
0,08	0,0012	0,008075541
0,16	0,00503475	0,0157619
0,23	0,0121	0,02953131
0,31	0,02278947	0,04097326
0,39	0,0349	0,05493775
0,47	0,05233831	0,07522751
0,55	0,0775	0,10523
0,63	0,1089077	0,149089
0,70	0,158	0,20752
0,78	0,2199285	0,2783817
0,86	0,2868	0,37033
0,94	0,3583839	0,469151

Н опт=0,07825423

Математична основа способу Мілна

Для реалізації засобу Мілна необхідно мати інформацію о попередніх точках. Тому засіб Мілна реалізуєтся після підрахунків по засобу Рунге-Кутта з заданной точностью.

Формула прогнозу:

Визначаємо значення проізводної:

Формула коррекції:

Якщо, то

Якщо, то закінчуємо розрахунок і даємо поманду на друк результата.

8.Блок-схема реалізації головного модуля, програми та графічної частини за засобом Мілна

нет

1.9.Реалізація програми за способом Мілна

DECLARE SUB MILN (T!, X!, Y!, A%, B%, H!, N%, E!, C%, X(), Y(), T(), LX3!, LY3!, LX2!, LY2!, LX1!, LY1!, XP!, YP!, XK!, YK!, MPX!, MPY!, MKX!, MKY!, XK1!, YK1!)

DECLARE SUB KUTT (T!, X!, Y!, A%, B%, H!, N%, E!, C%, X(), Y(), T(), KX1!, KY1!, KX2!, KY2!, KX3!, KY3!, KX4!, KY4!)

DECLARE SUB GRAF (T!, X!, Y!, A%, B%, H!, N%, E!, C%, X(), Y(), T(), LX3!, LY3!, LX2!, LY2!, LX1!, LY1!, XP!, YP!, XK!, YK!, MPX!, MPY!, MKX!, MKY!, XK1!, YK1!)

INPUT "C"; C%

E! = C% * 10 ^ (-4)

H! = E! ^ (1 / 4)

CONST A% = 0: CONST B% = 1

DIM SHARED T!(2000), X!(2000), Y!(2000), KX1!(2000), KY1!(2000), KX2!(2000), KY2!(2000), KX3!(2000), KY3!(2000), KX4!(2000), KY4!(2000)

DIM SHARED LX1!(2000), LY1!(2000), LX2!(2000), LY2!(2000), LX3!(2000), LY3!(2000), XP!(2000), YP!(2000), XK!(2000), YK!(2000), MPX!(2000), MPY!(2000), MKX!(2000), MKY!(2000), XK1!(2000), YK1!(2000)

T(0) = 0: X(0) = 0: Y(0) = 0

CALL KUTT(T!, X!, Y!, A%, B%, H!, N%, E!, C%, X(), Y(), T(), KX1!, KY1!, KX2!, KY2!, KX3!, KY3!, KX4!, KY4!)

FOR I% = 0 TO N%

PRINT T(I%), X(I%), Y(I)

NEXT I%

INPUT L!

CALL GRAF(T!, X!, Y!, A%, B%, H!, N%, E!, C%, X(), Y(), T(), LX3!, LY3!, LX2!, LY2!, LX1!, LY1!, XP!, YP!, XK!, YK!, MPX!, MPY!, MKX!, MKY!, XK1!, YK1!)

INPUT P!

CALL MILN(T!, X!, Y!, A%, B%, H!, N%, E!, C%, X(), Y(), T(), LX3!, LY3!, LX2!, LY2!, LX1!, LY1!, XP!, YP!, XK!, YK!, MPX!, MPY!, MKX!, MKY!, XK1!, YK1!)

FOR I% = 3 TO N%

IF (MPX(I% + 1) - MKX(I% + 1)) > E! THEN

XK1(I% + 1) = X(I% - 1) + (1 / 3) * H! * (MKX(I% + 1) + 4 * LX3(I%) + LX2(I% - 1))

PRINT T(I% + 1), XK1(I% + 1)

ELSE

XP(I% + 1) = XK(I% + 1)

PRINT T(I% + 1), XK(I% + 1)

END IF

FOR I% = 3 TO N%

IF (MPY(I% + 1) - MKY(I% + 1)) > E! THEN

YK1(I% + 1) = Y(I% - 1) + (1 / 3) * H! * (MKY(I% + 1) + 4 * LY3(I%) + LY2(I% - 1))

PRINT T(I% + 1), YK1(I% + 1)

ELSE

YP(I% + 1) = YK(I% + 1)

PRINT T(I% + 1), YK(I% + 1)

END IF

NEXT I%

INPUT M!

CALL GRAF(T!, X!, Y!, A%, B%, H!, N%, E!, C%, X(), Y(), T(), LX3!, LY3!, LX2!, LY2!, LX1!, LY1!, XP!, YP!, XK!, YK!, MPX!, MPY!, MKX!, MKY!, XK1!, YK1!)

END

SUB GRAF (T!, X!, Y!, A%, B%, H!, N%, E!, C%, X(), Y(), T(), LX3!, LY3!, LX2!, LY2!, LX1!, LY1!, XP!, YP!, XK!, YK!, MPX!, MPY!, MKX!, MKY!, XK1!, YK1!)

SCREEN 2

VIEW (170, 50)-(470, 150)

WINDOW (-1, 1)-(1, -1)

FOR I% = 0 TO N% - 1

PSET (T(I%), KX(I%))

PSET (T(I%), YK(I%))

LINE (T(I%), XK(I%))-(T(I% + 1), XK(I% + 1))

LINE (T(I%), YK(I%))-(T(I% + 1), YK(I% + 1))

NEXT I%

FOR I% = 4 TO N% - 1

PSET (T(I%), KX(I%))

PSET (T(I%), YK(I%))

LINE (T(I%), XK(I%))-(T(I% + 1), XK(I% + 1))

LINE (T(I%), YK(I%))-(T(I% + 1), YK(I% + 1))

PSET (T(I%), XK1(I%))

PSET (T(I%), YK1(I%))

LINE (T(I%), XK1(I%))-(T(I% + 1), XK1(I% + 1))

LINE (T(I%), YK1(I%))-(T(I% + 1), YK1(I% + 1))

NEXT I%

LINE (-1, 0)-(1, 0)

LINE (0, -1.5)-(0, 1.5)

END SUB

SUB KUTT (T!, X!, Y!, A%, B%, H!, N%, E!, C%, X(), Y(), T(), KX1!, KY1!, KX2!, KY2!, KX3!, KY3!, KX4!, KY4!)

N% = (B% - A%) / H!

I% = 0

T(I%) = T(0) + I% * H!

KX1(I%) = H! * (-2 + X(I%) + 5 * Y(I%))

KY1(I%) = H! * ((EXP(.5 * Y(I%) + T(I%)) - EXP(-.5 * Y(I%) + T(I%))) / 3 +.5 * Y(I%))

KX2(I%) = H! * (-2 + (X(I%) + KX1(I%) / 2) + 5 * (Y(I%) + KX1(I%) / 2))

KY2(I%) = H! * ((EXP(.5 * (Y(I%) + KY1(I%) / 2) + (T(I%) + H! / 2) - EXP(-.5 * Y(I%) + KY1(I%) / 2 - (T(I%) + H! / 2))) / 3 +.5 * (Y(I%) + KY1(I%) / 2)))

KX3(I%) = H! * (-2 + (X(I%) + KX2(I%) / 2) + 5 * (Y(I%) + KY2(I%) / 2))

KY3(I%) = H! * ((EXP(.5 * (Y(I%) + KY2(I%) / 2) + (T(I%) + H! / 2) - EXP(-.5 * Y(I%) + KY2(I%) / 2 - (T(I%) + H! / 2))) / 3 +.5 * (Y(I%) + KY2(I%) / 2)))

KX4(I%) = H! * (-2 + (X(I%) + KX3(I%) / 2) + 5 * (Y(I%) + KY2(I%) / 2))

KY4(I%) = H! * ((EXP(.5 * (Y(I%) + KY3(I%) / 2) + (T(I%) + H! / 2) - EXP(-.5 * Y(I%) + KY3(I%) / 2 - (T(I%) + H! / 2))) / 3 +.5 * (Y(I%) + KY3(I%) / 2)))

X(I% + 1) = X(I%) + (1 / 6) * (KX1(I%) + 2 * KX2(I%) + 2 * KX3(I%) + KX4(I%))

Y(I% + 1) = Y(I%) + (1 / 6) * (KY1(I%) + 2 * KY2(I%) + 2 * KY3(I%) + KY4(I%))

I% = I% + 1

LOOP UNTIL I% > N%

END SUB

SUB MILN (T!, X!, Y!, A%, B%, H!, N%, E!, C%, X(), Y(), T(), LX3!, LY3!, LX2!, LY2!, LX1!, LY1!, XP!, YP!, XK!, YK!, MPX!, MPY!, MKX!, MKY!, XK1!, YK1!)

N% = (B% - A%) / H!

FOR I% = 3 TO N%

T(I%) = T(0) + I% * H!

LX3(I%) = -2 * X(I%) + 5 * Y(I%)

LY3(I) = (EXP(.5 * (Y(I%) + T(I%)) - EXP(-.5 * Y(I%) - T(I%))) / 3 +.5 * Y(I%))

LX2(I% - 1) = -2 * X(I% - 1) + 5 * Y(I% - 1)

LY2(I% - 1) = (EXP(.5 * (Y(I% - 1) + T(I% - 1)) - EXP(-.5 * Y(I% - 1) - T(I% - 1))) / 3 +.5 * Y(I% - 1))

LX1(I% - 2) = -2 * X(I% - 2) + 5 * Y(I% - 2)

LY1(I% - 2) = (EXP(.5 * (Y(I% - 2) + T(I% - 2)) - EXP(-.5 * Y(I% - 2) - T(I% - 2))) / 3 +.5 * Y(I% - 2))

XP(I% + 1) = X(I% - 3) + (4 / 3) * H! * (2 * LX3(I%) - LX2(I% - 1) + 2 * LX1(I% - 2))

YP(I% + 1) = Y(I% - 3) + (4 / 3) * H! * (2 * LY3(I%) - LY2(I% - 1) + 2 * LY1(I% - 2))

MPX(I% + 1) = 2 + XP(I% + 1) + 5 * YP(I% + 1)

MPY(I% + 1) = EXP(.5 * YP(I% + 1) + T(I% + 1)) - EXP(-.5 * YP(I% + 1) - T(I% + 1)) / 3 +.5 * YP(I% + 1)

XK(I% + 1) = X(I% - 1) + (1 / 3) * H! * (MPX(I% + 1) - 4 * LX3(I%) + LX2(I% - 1))

YK(I% + 1) = Y(I% - 1) + (1 / 3) * H! * (MPY(I% + 1) - LY3(I%) + LY2(I% - 1))

MKX(I% + 1) = -2 + XK(I% + 1) + 5 * YK(I% + 1)

MKY(I% + 1) = (EXP(.5 * YK(I% + 1) + T(I% + 1)) - EPX(-.5 * Y(I + 1) - T(I% + 1)) / 3 +.5 * YK(I% + 1))

NEXT I%

END SUB

1.10.Результати реалізації програми за способом Мілна

Т	Х	У
0	0	0
0,1565085	0,002034751	0,008549757
0,3130169	0,01578947	0,03372118
0,4695254	0,05033831	0,07872751
0,6560338	0,1149077	0,1481089
0,7825423	0,2199285	0,2483817
0,9390508	0,3783839	0,389151