|
Главная -> Справочник по алгоритмам 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,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 |
|