Доставка цветов в Севастополе: SevCvety.ru
Главная -> Справочник по алгоритмам

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 [36] 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78

le PRINTРЕШЕНИЕ СИСТЕМЫ ДИФФЕРЕНЦИАЛЬНЫХ

20 PRINTУРАВНЕНИЙ МЕТОДОМ РУНГЕ-КУТТА

30 inputBBEAHTE ЧИСЛО УРАВНЕНИЙ N=N

40 DIM VaO/ACNbKCN/fcns/WaO

50 ШРиТЗАДАйТЕ ШАГ Н=Н

60 inputЗАДАЙТЕ НАЧАЛЬНОЕ Х0=К

70 FOR J=l ТО N!print!2.0(ЗАДАЙТЕ НАЧАЛЬНОЕ VOCJ)

80 input wc.osletv<:j)=w<:j>:next j

90 БОВиВ 200:FOR J=l TO N!LETU=H)ttF<:J) 100 LETKCJ>=U;LETV<:J>=W<:.J)+U/2:NEXT j 110 LETX=X+H/2:60SUB 260!FOR J=l TO N 120 LETU=H»F<:J):LETK<J)=K<J)+2»U 130 LETV<:J>=UKJ)+U/2:NEXT j

140 60sub 288!f0r j=l to n!letu=h»f<:j) 158 letk<:j)=k<:j)+2»u:letv<:j>=w<:j>+U!NEXT j

168 LETX=x+H/2!pRINT!Fi.9!An3 X=y:GOSUB 288

178 FOR J=l TO N:letv<:J)=W<:j> + <:k<:J>+Hii!F<:j)>.S

188 PR1HT!2.8!v<:J>= !F1.9!V<:J>

198 leti.j<:.-0=v<j>!next JsGOTO 98

288 LETFa>=v<:2)!letf<:2) = <:v<l)/x-v<:2>)/x-va>

218 return:END.

Пример. Для системы уравнений (см. подпрограмму, записанную со строки 200)

вида

при Л/ = 2, А = 0,1, хс = 0,2, 4/10 = 0,099500833 и у2с = 0,49235 будем получать следующие значения х, у, и ys

X у,

0,3 0,1483051278

0,4 0,1960019791

0,5 0,2422342383

1 0,4399784882

1,5 0,5578420795

0,4831085139 0,4702289086 0,453843868

0,3250861127

0,1398423554

y"=~=F(x,y,y), .

имеющий погрешность R-ih), реализуется с помощью следующих формул [36]: K,=hF (хг,у,;у,)\

К. = кр{х+:у,+\у<+Кй ;

К,=Нр{х1+\,,у1+\у.+\кйУ.+)-, К,=hF (х,+Л; у1+ hy!+а /(з; t + Кз) ;

i/i+, =«/, + /! у;+1-(к + К2 + Кз)

0,5766255981

-0,06446334779

Ввиду особого значения и широкого применения дифференциальных уравнений второго порядка полезно иметь специальную программу для их решения.

Метод Рунге - Кутта для дифференциального уравнения второго порядка

i/f-n =i/f+-g-(1+2K2-f 2Кз + /(4).

Перед началом вычислений надо задать шаг k и начальные значения хо, у (хо) ={/п и у (хс) =уо-

Программа 4.52.

18 printрешение дифференциального уравнения второго 28 print порядка методом рунге-кутта 38 inputbbeahte шаг н=н!inputвведите начальное х0=х 40 inputведите начальное у0=ы

50 inputbbeahte начальное dv/dx0=u!letv=w!LETZ=u 60 60sub 150:LETA=H»F:LETX=X+H/2 70 иЕТУ=ы+и*н/2+АжН/8! letz=u+a/2 80 60sub 150:LETB=H»F:LETZ=U+B/2

90 60sub 158!letc=h!kf!letx=x+h/2 100 letv=w+h»u+h3i!c/2!letz=u+c! gosub 158 118 letv=u+h»<:u+<:a+b+c>/6>!LEtw=v 120 letz=u+<a+<b+c)*2+h*f>6!letu=z 130 рр1итдля x=x!printv=v

140 printdv/r«=z:goto 68

158 letf=-v+<l-v*v)*z)tt28!return:end

11.4



Пример. Для дифференциального уравнения (см. подпрограмму, записанную в строке 150)

£={/"=-i/-f {20у)

при /i=s0,003125, Хо=0. {/0 = 2,0077 н 1/6 = 0 будем получать:

X у у-

0,003125 2,007690786 - 0,00571555801

0,00625 2,007665302 - 0,01044479113

6,009375 . 2,007626354 - 0,01435793978

0,0125 2,00757266 - 0,01759589125

0,015625 2,007516959 - 0,02027523125

Автоматическое изменение шага в ходе решения систем дифференциальных уравнений необходимо, если решение требуется получить с заданной точностью. При высокой точности (погрешность е=£=ЫО-) и решении в виде кривых с сильно различающейся крутизной автоматическое изменение шага обеспечивает уменьшение общего числа шагов в несколько раз, резко уменьшает вероятность возникновения числовой неустойчивости, дает более равномерное расположение точек графика кривых (решений) при их выводе на печать.

Метод Рунге - Кутта с автоматическим изменением шага заключается в том, что после вычисления )/(,+1) с шагом h все вычисления проводятся повторно с шагом Л/2. Полученный результат Jfa+i) сравнивается с Yi(i+,y Если I К*(,+ , <£, вычис-

ления продолжают с шагом h, в противном случае шаг уменьшают. Если это неравенство слишком сильное, шаг," напротив, увеличивают. При той же погрешности ~</г) лучшие результаты дает описанный ниже метод.

Метод Рунге - Кутта - Мерсона с автоматическим изменением шага обеспечивает приближенную оценку погрешности на каждом шаге интегрирования. Погрешность интегрирования имеет порядок . А° [3.1]. Этот метод реализуется следующим а[лго-ритмом.

1. .Задаем число уравнений N, погрешность 8 = £, начальный шаг интегрирования /г = Я и начальное значение хо, yio, Ут-

2. С помощью 5 циклов с управляющей

переменной У=1, 2.....Л/ вычисляем ко:*ф-

фициенты:

Ko, = AF, (X,; Yi,):

A:,i=/iF;(x,+-A;K(,+- ; Л..-.. Кт, = hFi (х,- + h: Y,i + -g- Ко Си) ; K3,=hFi {xi+~h- Yii+-L Кщ+ K2): K,ihFi(Xi + h: y,.,4-i-Ko,---/C2,-t-2K:<,) .

3. Находим (в последнем цикле) значение

Yi(,+ ,)= Ya-h (Д:п;+4А:з, + К.;) /6 и погрешность

/?д,+1)=(-2Ко,+9д:2/-8д:з,+а:4,)/зо. .

4. Проверяем выполнение условий

1/?и,-+1)1£/30.

Если первое условие ие выполняется, де лим шаг /г на 2 и повторяем вычисления с п. 2, восстановив .нанальные значения К,*. Если это условие выполняется и выполняется второе условие, значения x,+ i=Xi4-/i и ),(, +о выводятся на печать. Если второе условие ие выполняется, шаг h увеличивается вдвое и вычисления опять . повторяются с п. 2.

Таким образом, K,(„+i) выводится на печать только при одновременном выполнении условий этого пункта.

Программа 4.53.

Пример. Для системы дифференциальных уравнений

/Г,. = .=«/,+у2 Х» + Х2,.

dyi dx

= -2i/,-f4j/2+2x=-4x-7

при Л/ = 2, fe = 0,l, £=1.10-, хо = 0, 4/10=0 и у20 = 2 получим:

0,025

6,25-10-"

2,025

0,05

2,5-10-

2,05

0,075

5,625-10-

2,075

1-10-=

Следовательно, в- данном случае имеет место автоматическое уменьшение шага. Для /г = = 0,02 и £=Ы0- получим:

X f/l </2

0,02 6,731133879-10-* 2,020963777

0,04 1,898781104 •10- 2,041032153

ЬМ 6,756088366-10- 2,081182964

0,16 2.669871963-10-= 2,161549932

0,32 1,033398344 • 10 2,322641762

0,64 4,127175414-10-- 2,74775972

0,96 9,309436052-10" 2,981657806

1,28 . 1,665020406 3,338953272

В данном случае шаг автоматически увеличивается от значения /г = 0,02до ./г>==0,32. Время вычислений на одном., шаге (кроме начального при. уменьшении h) составляет около 5 с. • • -

Подпрограмма вычисления производных записывается со строки 400.

Как отмечалось, погрешность Нц,+ \) на каждом шаге метода Рунге - Кутта - Мерсона оценивается прибли.женно. При решении нелинейных дифференциальных уравнений истиннаи погрешность может отличаться . в несколько раз от заданной £.



Метод Рунге - Кутта - Фельберга с автоматическим изменением шага дает более точную оценку погрешности на каждом шаге и реализуется последовательным циклическим вычислением по следующим формулам (31):

Koi = hFi(x/; Ksi = /!Fj(xi+i-/u y,+±-K„i+-K>i);

/С„=hF, (x,+Л; к, -i Кщ +~Км~

, = hF,(xi+h;Y,--K.,--K,+

Yi(H-i)=yii+ X<!+ Xi XIJ2

(4.37)

Погрешность

I 2 J g

Kij-

в этом методе - разность приращений Ук+О. вычисленных по двум формулам: порядка п = 4 (4.37) и порядка « + Последняя формула не приводится, но использована для вычисления Рщ+п. Если У?у(,ч-1)> £. шаг h уменьшается вдвое, если Ri,i + ,<:E/20, он увеличивается вдвое. Этот метод имеег четвертый порядок.

05 PRINT-РЕШЕНИЕ СИСТЕМЫ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ 10 PRINTМЕТОДОМ РУНГЕ-КУТТА-МЕРСОНА С АВТОМАТИЧЕСКИМ 15 PRINT ВЫБОРОМ ШАГА:PRint!2.0!

48 INPUTBBEAHTE ЧИСЛО УРАВНЕНИЙ N=N

45 ИМ уао/ыао..А<:н)гС<н)/г1<:н)/Е<н>,Еао

50 inputЗАДАЙТЕ ПОГРЕШНОСТЬ ВЫЧИСЛЕНИЙ Е=Е1

68 inputЗАДАЙТЕ НАЧАЛЬНЫЙ ШАГ Н=Н

?0 INPUTЗАДАЙТЕ НАЧАЛЬНОЕ Х0=Х

80 FOR J=l TQ NSPRINTВВЕДИТЕ НАЧАЛЬНОЕ V0<J>=

90 INPUT WCJJrsLETVCJJWCJJrsNEXT j

110 LETE3=8s GOSUB 400sLETri=0sFOR J=l TO N

120 leta<:j>=F<:j))ttHsLETV<:j)=W<J)+a<:j)/3!next J

130 LETX=X+H/3!60SUB 4005FOR J=l TO N

140 LETV<J>=W<J>+(A<:J>+F<J>3I!H)-6sNEXT j

150 GOSUB 400SFOR J=l TO NsLETC<J)=F<J>»H

160 LETV<:J)=W<J>+A<J>/8+.375»C<:J>!№XT j

178 LETX=X+H.6s GOSUB 408: FOR J=l TO N

180 letd<:j)=f<:j)»H!LETV<:j)=w<j)+a(J)/2-i.5«c<:.J)+2«d<..o

190 next j!letx=X+H-2! GOSUB 408

200 FOR j=i TO H5LETe<:j)=f<:j)»h

210 letv<j>=w<J)+<:a<:j)+4»d<j>+e<:j>)/6

228 LETE£=ABS<:-2iifA<:J)+9»c<:J)-S)ttri<:J)+E<:J))30

230 IF E2<=E1 then 250

240 LETE3=lsG0T0 260

256 IF E2<El/20 THEN LETII=ri+l

260 next JsiF E3=0 then 296

276 LErx=X-HsFOR .J=l TO NsLETVCJ>=W«:JjsNEXT J

280 leth=H/2sG0T0 118

290 IF ri=N then leth=H+H

3C10 РР1НТ!Е1.9!ДЛЯ X=X5F0R J=l TO N

318 PRINT!2.0!V<:J)= !Fi.9!v<:J)

320 letw<:j>=v<:j>!Next jsgoto 110 400 iETFa>=va>+v<:2)-x»x+x-2

410 letf<:2>=-2»va)+4)ttV<:2>+2«X»X-4»X-7 420 RETLlRNsEND



0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 [36] 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78



0.4382
Яндекс.Метрика