численные методы. лекции чм. 1. Постановка задачи Коши
 Скачать 4.45 Mb.
|
|
8.Вложенный метод Рунге-Кутта При численном решении задачи Коши (1.1) - (1.2) иногда возникает проблема выбора длины шага. В этом случае часто применяют вычислительные схемы, основанные на нескольких методах Рунге-Кутта, что позволяет создавать компьютерные программы с автоматическим выбором длины шага. Фельберг (1969) предложил применить метод Рунге-Кутта пятого порядка точности ( р = 5 ), для оценки погрешности метода четвертого порядка точности ( р = 4 )таким образом, что на каждом шаге требуется шесть раз вычислять значение f (в то время, как применение формулы Рунге-Ромберга (5.4) для классического метода Рунге-Кутта требует вычисления f в восьми точках). Соответствующий вложенный метод Рунге-Кутта определяется формулами : ( 8.1 ) ( 8.2 ) где ¢ = ¢ = ¢ = ¢ = - - ì í ï ïï î ï ï ï y y y y y y y x y x y 1 2 2 3 3 4 4 4 2 3 4 2 y x y x y x y x 1 2 3 4 1 1 0 1 1 2 = = = = = = = = , X K Y K y x K ( ) ( ) x y 2 0 ¢¢ ¢¢ = e y x y k k k = - £ - ( ) 10 5 y y h c K c K k k + = + + + 1 1 1 6 6 ( ) y y h c K c K k k * * * ( ) + = + + + 1 1 1 6 6 22 ( 8.3 ) Коэффициенты выбираются так, что формулы (8.1) и (8.2) задают методы пятого и четвертого порядков точности соответственно. При использовании формул ( 8.1 ) - ( 8.3 ) необходимо учесть , что для получения на каждом последующем шаге решения более низкого порядка точности - , используется решение на предыдущем шаге, имеющее более высокий порядок точности. Кэш и Кэрп (1992) показали, что ошибки округления минимизируются, если значения коэффициентов берутся согласно таблицы 7. Пусть - решение задачи Коши, вычисляемое по формуле (8.1) , - решение, вычисляемое с помощью формулы ( 8.2 ). Для контроля погрешности, на каждом шаге необходимо вычислять: ( 8.4 ) Если положить ,где - точное решение задачи Коши, - решение по формуле (8.2),тогда можем записать, что . Тогда с учетом формулы (2.6) можем записать , где р = 5 ( 8.5 ) Таблица 7 i 1 2 1/5 1/5 0 0 3 3/10 3/40 9/40 4 3/5 3/10 -9/10 6/5 K f x y K f x a h y b K K f x a h y b K b K k k k k k k 1 2 2 21 1 6 6 61 1 65 5 = = + + = + + + + ( , ) , ( , ) , . . . . . . . . . . . . . . . . . . . . . . . . . . , ( , ) a b c c i i j i i , , , , * y k + 1 * y k y x k k ( ) y x k k * ( ) D( ) ( ) ( ) ( ) * * h y h y h c c K k k i i i i = - = - + + = å 1 1 1 6 e y x y k k k * * ( ) = - y x k ( ) y k * e h k * ( ) » D D( ) h Ch p » a i j b i j c i c i * 37 378 2825 27648 250 621 18575 48384 125 594 13525 55296 23 5 1 -11/54 5/2 -70/27 35/27 0 6 7/8 1/4 В таблице 7 индекс j изменяется в пределах Предположим, что - значение шага h, при котором достигается решение данной задачи Коши с требуемой точностью . Положим , тогда с учетом ( 8.5 ) . Для любого другого значения соответственно запишем . Тогда из отношения погрешностей получим ( 8.6 ) Конечно, вначале значения не известны. На самом деле, нет необходимости задавать , если известны и . Поэтому значение приходится выбирать эмпирически. Например, можно связать величины и получив численное решение одним шагом метода Эйлера . Тогда, сделав шаг мы можем оценить, насколько численное решение с этим шагом - удовлетворяет требуемой точности. Так, если то отношение ( 8.6 ) показывает, во сколько раз нужно пересчитать (уменьшить) неудачный выбранный шаг . С другой стороны, если отношение ( 8.6 ) покажет нам, во сколько раз мы можем безопасно увеличить удачно выбранный шаг Если выполнено условие то для следующего пробного значения h= можем записать ( 8.7 ) Аналогично, получаем выражение и для других пробных шагов: и т.д. , j=1,…N ( 8.8 ) При использовании формул ( 8.7 ) - ( 8.8 ) необходимо учитывать следующие особенности. Так как зависит от шага h, то уменьшение шага ( в случае ) говорит о уменьшении порядка точности оценки , т.е. показатель 277 14336 1631 55296 175 512 575 13824 44275 110592 253 4096 512 1771 1 1 £ £ - j i h 0 e D D 0 0 = ( ) h D 0 0 5 » Ch h 1 D D D 1 1 1 1 5 = » ( ) , h Ch h h 0 1 0 1 0 2 » æ è ç ö ø ÷ D D h 0 0 , D h 0 D 0 h 1 D 0 D 0 h 1 { } D 0 0 1 0 0 » + e y h f x y ( , ) h 1 y 1 D D 1 0 > h 1 D D 1 0 < h 1 D D 1 0 < h 2 h h 2 1 0 1 0 2 » æ è ç ç ö ø ÷ ÷ D D h h 3 2 0 2 0 2 » æ è ç ç ö ø ÷ ÷ D D h h j j j » - - 1 0 1 0 2 D D D 0 D D 1 0 > D 1 24 степени в формуле ( 8.5 )заменяется на р=4. Соответственно показатель степени в формуле ( 8.6 ) 0.2 заменяется на 0.25 и следующий шаг рассчитывается исходя из этого обстоятельства. Иными словами, можно записать: ( 8.9 ) На практике для повышения эффективности вычислений вводят некоторый множитель “безопасности” , который на несколько процентов меньше единицы - S 0.9. ( 8.10 ) Коротко рассмотрим алгоритм метода для его компьютерной реализации. Обозначим численное решение получаемое методом Рунге-Кутта с переменным шагом (j=1,2,...N) в не равноотстоящих узловых точках через . Введем также параметр , характеризующий точность вычислений - . Для методов Рунге- Кутта четвертого порядка точности достаточно принять Погрешность на каждом шаге определим следующим образом: , k=0,...,N-1 ( 8.11 ) 1. Задаются начальные параметры , ,пробный шаг , по формуле (8.11) вычисляется погрешность = . Для заданного пробного шага h= в точке = по формуле (8.3) последовательно вычисляются коэффициенты ( I=1,..6 ),а по формуле (8.4 )положив h= рассчитывается ошибка Если выполнено условие | / |<1 , то шаг выбран удачно. Если же окажется, что | / |>1 ,то шаг считается неудачным и пересчитывается по формуле ( 8.10 ) h h h 0 1 0 1 0 2 1 0 1 0 1 0 25 1 0 = £ > ì í ï ï î ï ï D D D D D D D D » h Sh Sh 0 1 0 1 0 2 1 0 1 0 1 0 25 1 0 = £ > ì í ï ï î ï ï D D D D D D D D h j x x h j j j = + -1 y j e e £ - 10 4 D 0 { } D 0 ( ) ( , ) h y h f x y j k j k k » + e x 0 y x ( ) 0 h 1 D 0 { } e y x h f x y x ( ) ( , ( ) 0 1 0 0 + h 1 x 1 x h 0 1 + K i h 1 D 1 D 1 D 0 ( ) D D 1 0 < D 1 D 0 ( ) D D 1 0 > h 1 25 С новым значением шага , вычисления повторяют, получая новое значение Далее по формуле ( 8.10 ) вычисляется значение следующего пробного шага: и по формуле ( 8.2 ) определяется численное решение задачи Коши в точке Таким образом, на этом этапе определяется численное решение с шагом в точке ,а также прогнозируется следующее значение шага интегрирования - 2.Значение следующего шага h, принимают равным ,и по формуле(8.11) подсчитывают величину = а в точке рассчитывается ошибка Если шаг неудачный, ,он пересчитывается: Значение следующего пробного шага полагают равным На этом этапе определяется численное решение в точке , а также прогнозируется следующее значение шага интегрирования - 3.Следующий шаг принимается равным и в точке рассчитывается ошибка и , а описанная выше процедура повторяется. На этом этапе получают численное решение в точке Таким образом проходим весь интервал интегрирования [a, b] и получаем набор точек и численное решение в них - (j=1,2,…,N). Особое внимание следует обратить на то, чтобы при решении задачи Коши избежать перехода через конечную точку интервала интегрирования - . Для этого в процессе вычислений необходимо следить за знаком следующей величины: ! h sh 1 1 1 0 0 25 = - D D h h 1 1 = ! D 1 h sh 2 1 1 0 0 2 = - D D x 1 y 1 h 1 x 1 h 2 h 2 D 0 ( ) { } e y h f x y 1 2 1 1 + , x x h 2 1 2 = + D 2 D D 2 0 1 > ! h sh 2 2 2 0 0 25 = - D D h sh 3 2 2 0 0 2 = - D D y 2 x 2 h 3 h h = 3 x x h 3 2 3 = + D 0 D 3 y 3 x 3 x x h j j j = + -1 y j X N 26 , (j=1,...,N) Если t<0 , то текущая точка интегрирования находится внутри интервала интегрирования , если t>0 то текущая точка интегрирования находится вне этого интервала. В случае t>0 можем записать, что , тогда, положив, по формуле (8.2) получим численное решение с шагом в конечной точке интегрирования При решении задачи Коши для системы из m дифференциальных уравнений формулы (8.4) и (8.11) соответственно записываются следующим образом ; где i=1,..,m , k=0,..,N-1 , а коэффициенты вычисляются для системы из m уравнений. Формула (8.8) записывается так C учетом этого записывается и формула (8.10) Пример 6. Решением на отрезке [0,1] задачи Коши: ( 8.14 ) у(0) = 1 является функция Влияние члена содержащего экспоненту наиболее заметно на интервале ( 0, 0.2 ), для значений х > 0.2 его величина практически не сказывается на решении, и оно становится медленно меняющимся. Теперь рассмотрим механизм управления длиной шага в методе Рунге-Кутта для этой задачи .Задав произвольно начальный шаг , определим исходя из его значения действительный шаг и t x h x x h x j j N j j = + - + - - - ( ) ( ) 1 1 0 [ ] x x N 0 , x h x N N N - + > 1 h x x N N N = - -1 y N h N X N D i i k i k n n n n y h y h c c K = - = - + + = å 1 1 1 6 ( ) ( ) ( ) * * ( ) D i i k j k k k mk y h f x y y y 0 1 2 = + e ( , , , . . , ) K n h h j j i m i i j = ì í ï îï ü ý ï þï æ è ç ç ö ø ÷ ÷ - £ £ - 1 1 0 1 0 2 max D D h Sh Sh i m i i i m i i i m i i i m i i 0 1 1 0 1 0 2 1 0 1 1 1 0 1 0 25 1 0 1 1 1 = × ì í î ü ý þ æ è çç ö ø ÷÷ ì í î ü ý þ > × ì í î ü ý þ æ è çç ö ø ÷÷ ì í î ü ý þ < ì í ï ïï î ï ï ï £ £ £ £ £ £ £ £ max , max max , max D D D D D D D D ¢ = - + + y x y x x x ( ) ( ) cos sin 25 25 y x x e x ( ) sin = + - 25 h 1 !h 1 27 вычислим последующий шаг . Вычисления будем проводить с различными значениями точности . Таблица 8 точность h1 h1^ h2 .3679E+00 .10 .100000 .500000 .1353E+00 .10 .100000 .500000 .4979E-01 .10 .100000 .500000 .1832E-01 .10 .100000 .500000 .6738E-02 .10 .100000 .500000 .2479E-02 .10 .100000 .447400 .9119E-03 .10 .100000 .366300 .3355E-03 .10 .100000 .299901 .1234E-03 .10 .100000 .245538 .4540E-04 .10 .100000 .201030 .1670E-04 .10 .100000 .164589 .6144E-05 .10 .100000 .134754 .2260E-05 .10 .100000 .110327 .8315E-06 .10 .100000 .090328 .3059E-06 .10 .070412 .082856 Как видно из таблицы 8 задание различных значений не слишком значимо и пробный шаг оказывается вполне удачным (значения пробного и действительного шага совпадают), что позволяет проводить вычисления уже с большим значением шага ( ). Однако по мере увеличения требований точности уменьшается и длина шага, так при принятую величину шага = 0.1 уже приходится пересчитывать. Решим задачу Коши (8.14) на всем интервале интегрирования методом Рунге-Кутта с переменным шагом. Примем требуемую точность вычислений = , начальную длину шага = 0.1.Результаты решения показаны на рис.9. Как видно из графика, численное решение методом Рунге-Кутта с переменной длиной шага полностью соответствует точному решению, в рамках принятой точности , при этом h 2 e e h 1 h 1 !h 1 h 2 e » - 3 10 7 * h 1 e 10 6 - h 1 e e k » - 10 6 28 Рис.7 Теперь решим задачу Коши методом Эйлера, с шагом h = 0.1 на отрезке [0,1]. Результаты решения показаны на рис.8. Из графика видно, что несоответствие численного и точного решений наблюдается на всем интервале интегрирования, а погрешность решения . Однако в случае уменьшения длины шага до h = 0.01 картина резко изменяется ( рис 9 ), при этом , что характерно для метода Эйлера. Графики представленные на рис.8 - рис.9 являются иллюстрацией того, как важно правильно выбрать длину шага при численном решении задачи Коши. Дальнейшее уменьшение длины шага h для метода Эйлера не приведет к уменьшению погрешности метода, а наоборот, может ухудшить численное решение ( см. рис 3.). e k >> 1 e k » × - 2 5 10 2 29 Рис.8 30 Рис.9 В заключении необходимо отметить, что дифференциальное уравнение в задаче ( 8.14 ) относится к так называемым жестким уравнениям . Жесткость - это свойство самого уравнения - в этом случае приходится вести численное решение с неоправданно мелким шагом, хотя само решение меняется медленно . Так, наличие малого параметра при старшей производной в дифференциальном уравнении, может служить признаком его жесткости. ПРИЛОЖЕНИЕ 1 Численные методы решения задачи Коши для обыкновенных дифференциальных уравнений реализованы во многих математических библиотеках . Одной из таких библиотек является математическая библиотека Numerical Recipes ( издательство Кембриджского университета )поставляемая вместе с компилятором Microsoft FORTRAN PowerStation. На примере этой библиотеки рассмотрим процедуру вызова подпрограммы решения задачи Коши (пример 6) методом Рунге-Кутта с переменным шагом (модуль xodeint), листинг которой приведен ниже. PROGRAM xodeint ****************************************************** * driver for Numerical Recipes Library. * 31 * Runge-Kutta method with stepsize control * ****************************************************** 1 INTEGER*4 KMAXX,NMAX,NVAR 2 INTEGER*1 scale, nx, ny 3 PARAMETER (KMAXX=200,NMAX=50,NVAR=1) 4 INTEGER i,kmax,kount,nbad,nok,nrhs 5 DOUBLE PRECISION dxsav,eps,h1,hmin,x1,x2,x,y,ystart(NVAR),testf,v 6 DOUBLE PRECISION TEST(KMAXX),RES(KMAXX) 7 CHARACTER string*30 8 COMMON /path/ kmax,kount,dxsav,x(KMAXX),y(NMAX,KMAXX) 9 COMMON nrhs 10 COMMON /param/ scale, string(3), nx,ny 11 EXTERNAL derivs,rkqs 12 testf(v)=sin(v)+exp(-25.*v) 13 open(8,file='hard1.res') 14 string(1) ='Stiff Equation' 15 string(2) ='Solutions' 16 string(3) ='Argument x' 17 scale = 9 18 nx = 4 19 ny = 4 20 nrhs=0 21** The integration cut *********************** 22 x1=0.0 23 x2=1.0 24** The boundary condition ******************** 25 ystart(1)= 1. 26********************************************** 27 eps=1.0D-6 28 h1= 0.01 29 hmin= 0.0 30 kmax=KMAXX !How many steps is stored 31 dxsav=0. !Interval for outputs 32 call odeint(ystart,NVAR,x1,x2,eps,h1,hmin,nok,nbad,derivs,rkqs) 33 do i =1, kount 34 35 test(i)= testf(x(i)) 36 res (i)= y(1,i) 37 enddo 38 write(*,'(/1x,a,t30,i3)')'Successful steps:',nok 39 write(*,'(1x,a,t30,i3)') 'Bad steps:',nbad 40 write(*,'(1x,a,t30,i3)') 'Function evaluations:',nrhs 41 write(*,'(1x,a,t30,i3)') 'Stored intermediate values:',kount 42 write(*,'(/1x,t9,a,t20,a,t33,a)') 'X','Y1','TEST' 43 write(8,'(/1x,a,t30,i3)') 'Successful steps:',nok 44 write(8,'(1x,a,t30,i3)') 'Bad steps:',nbad 45 write(8,'(1x,a,t30,i3)') 'Function evaluations:',nrhs 46 write(8,'(1x,a,t30,i3)') 'Stored intermediate values:',kount 47 write(8,'(/1x,t9,a,t20,a,t33,a)') 'X','Y1','TEST' 48 do i=1,kount 49 write(*,'(1x,f10.4,2x,2f14.6)') x(i),y(1,i),test(i) 50 write(8,'(1x,f10.4,2x,2f14.6)') x(i),y(1,i),test(i) 51 enddo 52 53 call grafik(x,res,test,kount) 54 close(8) 32 55 END 56*********************************************** 57* Right hand subroutine * 58*********************************************** 59 SUBROUTINE derivs(x,y,dydx) 60 INTEGER nrhs 61 DOUBLE PRECISION x,y(*),dydx(*) 62 COMMON nrhs 63 nrhs=nrhs+1 64 dydx(1)=-25.*y(1) + cos(x) + 25.*sin(x) 65 return 67 END В строке 3 задаются основные параметры NVAR - число уравнений в данной системе дифференциальных уравнений (NVAR=1), KMAXX-максимальное число узловых точек,по умолчанию принимается KMAXX=200, NMAX- максимальное число уравнений в системе, по умолчанию принимается NMAX=50.В строке 4 параметры kmax,kount,nbad,nok,nrhs определяют размерность вектора зависимой переменной (можно положить kmax=KMAXX), число точек вывода решения, число неудачных шагов, число удачных шагов, число обращений к подпрограмме правой части соответственно. В строке 5 параметры x1 и x2 задают начальную и конечную точку интегрирования соответственно. Параметр dxsav определяет шаг вывода решения dxsav=(x2-x1)/N, где N- число точек вывода. В зависимости от величины dxsav вывод решения действительно осуществляется в kount узловых точках, число которых kount . Если положить dxsav=0, то число точек и шаг вывода определяется самой программой. Параметры h1 и hmin задают начальный (пробный) шаг и минимальный шаг (обычно полагают hmin=0),массивы ystart и y определяют вектор начальных условий и матрицу вывода решения соответственно. Параметр testf определяет тестовое решение примера 6 через оператор-функцию в строке 12. Библиотека Numerical Recipes не содержит подпрограммы графического отображения информации, поэтому студенты могут написать ее сами, либо могут воспользоваться готовой подпрограммой grafik ,вызываемой в строке 53. В строке 6 массивы TEST и RES вводятся для построения графиков для тестирующего и численного решения. В строке 2 задаются параметры вывода графической информации - масштаб графика, сгущение сетки по оси х, сгущение сетки по оси у соответственно. Символьный массив string содержит название задачи, обозначение вертикальной и горизонтальной осей координат. Условие задачи Коши записывается с 21 по 25 строку, необходимая точность вычислений задается в строке 27. Вызов подпрограммы решения задачи Коши осуществляется в строке 32, подпрограмма записи правой части Derivs для данной задачи содержится в строках 56 - 67. £ N 33 В заключении покажем, как задавать исходные данные для решения примера 5 . Положим NVAR=4, ystart(1)=ystart(2)=0, ystart(3)=ystart(4)=2, x1=1.,x2=1.5, h1=0.1,hmin=0. Подпрограмма правой части SUBROUTINE derivs(x,y,dydx) INTEGER nrhs DOUBLE PRECISION x,y(*),dydx(*) COMMON nrhs nrhs=nrhs+1 dydx(1)=y(2) dydx(2)=y(3) dydx(3)=y(4) dydx(4)=-(4./x)*y(4)-(2./x/x)y(3) return END Библиотека Numerical Recipes собрана в файле nr32.lib, программа отображения графической информации находится в файле grsubs.for. Все данные с плавающей точкой необходимо задавать с двойной точностью. |