Главная Случайная страница


Категории:

ДомЗдоровьеЗоологияИнформатикаИскусствоИскусствоКомпьютерыКулинарияМаркетингМатематикаМедицинаМенеджментОбразованиеПедагогикаПитомцыПрограммированиеПроизводствоПромышленностьПсихологияРазноеРелигияСоциологияСпортСтатистикаТранспортФизикаФилософияФинансыХимияХоббиЭкологияЭкономикаЭлектроника






ЧИСЛЕННОЕ РЕШЕНИЕ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

Дифференциальное уравнение первого порядка,разрешенное относительно производной, имеет вид

y ' = f (x,y ). (6.1)

Решением дифференциального уравнения (6.1) называется функция φ (x), подстановка которой в уравнение обращает его в тождество: φ' (x)= f (x, φ (x)) . График решения y = φ ( x )называется интегральной кривой.

Задача Кошидля дифференциального уравнения (6.1) состоит в том, чтобы найти решение уравнения (6.1), удовлетворяющему начальному условию

(6.2)

Пару чисел (x0 , y0) называют начальными данными. Решение задачи Коши называется частным решениемуравнения (6.1) при условии (6.2).

Численное решение задачи Коши [формулы (6.1), (6.2)] состоит в том, чтобы получить искомое решение φ ( x )в виде таблицы его приближенных решений для заданных значений аргумента xна некотором отрезке [ a , b ]:

X0 = a , x 1 , x 2 , . . . , x m = b (6.3)

Точки по формуле (6.3) называются узловыми точками, а множество этих точек называется сеткойна отрезке [ a , b ] . Будем использовать равномерную сетку с шагом h :

h = ( b - a ) / m ; x i - x i - 1 = h или x i = x 0 + ih ( i = 1 , . . . , m ).

Приближенные значения численного решения задачи Коши в узловых точках x i обозначим через y i:

( i = 1 , . . . , m ).

Для любого численного метода решения задачи [см. формулы(6.1), (6.2)] начальное условие (6.2) выполняется точно , то есть .

Величина погрешности численного метода решения задачи Коши на сетке отрезка [ a , b ] оценивается величиной

.

Считается, что численный метод имеет pпорядок точности по шагу h на сетке , если расстояние d можно представить в виде степенной функции от h :

, p > 0,

где С- некоторая положительная постоянная, зависящая от правой части уравнения (6.1) и рассматриваемого метода. В данном случае очевидно, что когда шаг h стремится к нулю, погрешность d также стремится к нулю.

Метод Эйлера.Простейшим численным методом решения задачи Коши [см. формулы (6.1), (6.2)] является метод Эйлера. Угловой коэффициент касательной к интегральной кривой в точке P0 (x0 , y0) есть y '0 = f (x0 ,y0 ).

Найдем ординату y1 касательной , соответствующей абсциссе x 1 = x 0 + h . Так как уравнение касательной к кривой в точке P 0 имеет вид y - y0 = y ' ( x - x0 ) , то y1 = y0 + h f (x0 , y0 ).

Угловой коэффициент в точке P1 (x1 ,y1) также находится из данного дифференциального уравнения y'1 = f(x1,y1). На следующем шаге получаем новую точку P2 (x2 ,y2) , причем

x2 = x1 + h; y2 = y1 + hf(x1 ,y1) .

Продолжая вычисления в соответствии с намеченной схемой, получим формулы Эйлерадля m приближенных значений решения задачи Коши с начальными данными (x0 , y0) на сетке отрезка [ a , b ] с шагом h:

xi = xi -1 + h; yi = yi - 1 + hf(xi - 1 , yi - 1 ); где ( i = 1,2 , . . . , m ) (6.4)

Графической иллюстрацией приближенного решения (6.4) является ломаная, соединяющая последовательно точки P0 , P1, P2, . . . ,Pm , которую называют ломаной Эйлера .

Погрешность метода Эйлера можно оценить неравенством

 

или представить в виде

,где .

Это означает, что метод Эйлера имеет первый порядок точности. В частности, при уменьшении шага h в десять раз погрешность уменьшится примерно в десять раз.

Практическую оценку погрешности решения, найденного на сетке с шагом h/2, в точке производят с помощью приближенного равенства- (правила Рунге):

,(6.5)

 

где p- порядок точности численного метода. Таким образом, оценка полученного результата по формуле (6.5) вынуждает проводить вычисления дважды: один раз с шагом h , другой- с шагом h/2 .

Методы Рунге-Кутта.Численные методы решения задачи Коши

на равномерной сетке { x0 = a , x1 , x2 , . . . ,xm = b}отрезка[ a, b ]с шагом

h= (b -a)/m являются методами Рунге-Кутта, если, начиная с данных (x0 ,y0 ), решение ведется по следующим рекуррентным формулам:

(6.6)

Метод называют методом Рунге-Кутта порядка p, если он имеет p- й порядок точности по шагу h на сетке. Порядок точности p достигается с помощью формул (6.6) при определенных значениях коэффициентов cj и dj (j= 1,2, . . . ,p); c1 всегда полагают равным нулю.

Метод Рунге-Кутта второго порядка называют методом Эйлера-Коши, если p = 2 , c1 = 0 , c2 = 1 , d1 = d2 = 1/2. Алгоритм Эйлера-Коши получается из формул (6.6):

 

xi =xi-1 + h, yi = yi-1 + Δyi-1; Δyi-1 = (1/2)[ k1[i -1] + k2[i -1]] ( i = 1, ..., m); (6.7)

k1[ i - 1] = hf ( xi-1,yi-1), 2[ i - 1 ] = hf ( xi-1 + h , yi-1 + hf ( xi-1 ,yi-1) )

 

Для практической оценки погрешности решения можно применять правило Рунге, полагая в формуле (6.5) р = 2.

Метод Рунге - Кутта четвертого порядка называют классическим методом Рунге-Кутта, если p = 4, c1 = 0, c2 = c3 = 1/2, c4 = 1, d1 = d4 = 1/6, d2 = d3 =1/3.

Из рекуррентных формул (6.6)получим алгоритм решения задачи Коши классическим методом Рунге-Кутта :

 

x I = x i - 1 + h; y i = y i - 1 + Δy i – 1 ( i = 1,2, . . . , m);

 

Δyi-1 = 1/6 [ k1[ i - 1] + 2 k2[ i - 1] + 2k3[ i - 1] + k4[ i - 1] ];

 

k1[ i - 1] = h f ( xi - 1 , yi -1);

k2[ i - 1 ] = h f( xi - 1 + (1/2) h , y i - 1 + (1/2)k1[ i - 1 ] );

 

k3[ i - 1] = h f ( xi - 1 + (1/2)h , y i - 1 + (1/2)k2[ i - 1 ] );

 

k4[ i - 1 ] = h f( xi - 1 + h , y i - 1 + k3[ i - 1 ] ).

Графиком приближенного решения является ломаная, последовательно соединяющая точки Pi( xi , yi) ( i = 0 , 1 , 2, . . ., m). С увеличением порядка численного метода звенья ломаной приближаются к ломаной , образованной хордами интегральной кривой y = φ(x) , последовательно соединяющими точки ( xi , φ(xi) ) на интегральной кривой .

Правило Рунге (6.5) практической оценки погрешности решения для численного метода четвертого порядка

 

Задание:

Решить задачу Коши на равномерной сетке. Решение найти в четырех узловых точках (шаг h1 равен [ b - a ] / 4). Найти решение в тех же узлах, ведя расчет с уменьшенным вдвое шагом. Вычислить погрешности приближений при расчете с шагом h2 = h1 / 2.

Задачу решить с помощью системы MATHCAD следующим образом:

а) Эйлера;

б) Эйлера-Коши ;

в) Рунге-Кутта .

Варианты лабораторных работ:

1. ;

2. ;

3. ;

4. ;

5.

6. ;

7. ;

8. ;

9. ;

10. ;

11. ;

12. ;

13. ;

14. ;

15. ;

16. ;

17. ;

18. ;

19. ;

20. ;

Вспомогательные материалы

 

1. Пример:

Решить задачу Коши

; ; .

на равномерной сетке с шагом h = 0,1. Решение найти в четырех узловых точках.

С помощью программы найти решение в тех же узлах, ведя расчет с уменьшенным вдвое шагом. Вычислить погрешности приближений при расчете с шагом h = 0,05.

Задачу решить методами:

а) Эйлера;

б) Эйлера-Коши;

в) Рунге-Кутта.

Решение. Здесь f(x,y) = x + y ; m = 4 ; a = 0 ; b = 0,4 ;

h = ( b - a ) / m = 0,4 /4 = 0,1

 

1. Используя рекуррентные формулы x0 = 0 ; y0 = 1; xi = x i - 1 + 0,1; y i = y i - 1 + 0,1( x i - 1 + y i - 1 ) i = ( 1, 2, 3, 4 ) ,

последовательно находим

при i = 1 : x1 = 0,1 ; y1 = 1 + 0,1( 0 + 1 ) = 1,1;

при i = 2 ; x2 = 0,2 ; y2 = 1,1 + 0,1( 0,1 + 1,1 ) = 1,22;

при i = 3 ; x3 = 0,3 ; y3 = 1,22 + 0,1( 0,2 + 1,22 ) = 1,362;

при i = 4 ; x4 = 0,4 ; y4 = 1,362 + 0,1( 0,3 + 1,362 ) = 1,528 2.

С помощью программы находим решение при h = 0,05.

Обозначив , d i = | y i ( h ) - y i ( h/2) | сведем результаты вычислений в табл. 6.1

Таблица 6.1

I x i Y i ( h ) Y i ( h / 2) d
0,1 1,1 1,105 0,005
0,2 1,22 1,231 012 0,011 012
0,3 1,362 1,380 119 0,018 191
0,4 1,528 2 1,554 911 0,028 738

2. Формулы (6.7) в нашем случае принимают вид

 

k1[ i - 1] = h ( xi-1 + yi-1); k2[ i - 1 ] = h ( xi-1 + h + yi-1 + k1[ i - 1]);

xi =xi-1 + h; yi = yi-1 + (1/2)[ k1[i -1] + k2[i -1]] ( i = 1, 2 , 3 , 4).

 

Полагая x 0 = 0 , y 0 = 1 , последовательно находим

при i = 1:

 

k1[ 0 ] =0,1( 0 + 1 ) = 0,1; k2[ 0 ] = 0,1( 0 + 0,1 + 1 + 0,1) = 0,12;

x1= 0 + 0,1 = 0,1; y1 = 1 + ( 1/2 )( 0,1 + 0,12 ) = 1,11;

 

при i = 2:

k1[ 1 ] =0,1( 0,1 + 1,11 ) = 0,121; k2[ 1 ] = 0,1( 0,1+0,1+1,11+0,121) = 0,143 1;

x1= 0,1 + 0,1 = 0,2; y1 = 1,11+(1/2)(0,121+0,143) = 1,242 5.

Далее получаем при i = 3: x 3 = 0,3; y 3 = 1,398 465;

При i = 4: x 4 = 0,4; y 4 = 1,581 804.

С помощью программы проводим вычисления с половинным шагом. Результаты заносим в таблицу, аналогичную табл. 6.1.

3. Из формул (6.8) получаем

 

k1[ i - 1] = h ( xi-1 + yi-1);

k2[ i - 1 ] = h ( xi-1 + (1/2)h + yi-1 +(1/2) k1[ i - 1]);

k3[ i - 1 ] = h ( xi-1 + (1/2)h + yi-1 +(1/2) k2[ i - 1]);

k4[ i - 1 ] = h ( xi-1 + h + yi-1 + k3[ i - 1]);

xi =xi-1 + h , yi = yi-1 + (1/6)[ k1[i -1] + 2k2[i -1] + 2k3[i -1] + k4[ i - 1 ]]

для i = 1, 2 , 3 , 4.

Полагая x 0 = 0, y 0 = 1, последовательно находим

при i = 1:

k1[ 0 ] = 0,1(0 + 1 ) = 0,1;

k2[ 0 ] = 0,1(0 + 0,05 + 1 + 0,05) = 0,11;

k3[ 0 ] = 0,1(0 + 0,05 + 1 + 0,055) = 0,110 5;

k4[ 0 ] = 0,1(0 + 0,1 + 1 + 0,110 5) = 0,121 050;

x1= 0 + 0,1 = 0,1;

y1 = 1 +( 1/6 )( 0,1 + 2*0,11+2*0,110 5+0,121 05 ) = 1,110 342;

при i = 2 :

k1[ 1 ] = 0,1(0,1 + 1,110 342 ) = 0,121 034 2;

k2[ 1 ] = 0,1(0,1 + 0,05 + 1,110 342 + 0,060 517 1) = 0,132 638 5;

k3[ 1 ] = 0,1(0,1 + 0,05 + 1,110 342 + 0,066 042 95) = 0,132 638 5;

k4[ 1 ] = 0,1(0,1 + 0,1 + 1,110 342 + 0,132 638 5) = 0,144 298 0.

x2= 0,1 + 0,1 = 0,2;

y2 = y1 +( 1/6 ) [ k1[1] + 2 k2[1] +2 k3[1] + k4[1]] = 1,242 805.

Далее получаем при i = 3 x3 = 0,3; y3 = 1,399 717;

i = 4 x4 = 0,4; y4 = 1,583 648;

С помощью программы проводим вычисления с половинным шагом. Результаты заносим в таблицу.

Решение:

2. Блок-схема численного решения задачи Коши для дифференциального уравнения первого порядка методами Эйлера , Эйлера-Коши и Рунге-Кутта (рис 6.1)

Методом Эйлера
Методом Эйлера -Коши
    Методом Рунге -Кутта
Вычисление приближенного решения на одном шаге

Рис. 6.1

3. Пример программы для функции y / = x + y. (пример приведен для удобства разработки программы на MATHCADе студентами, привыкшими работать в ПАСКАЛЕ)).

program DifEquationsFirstOrder;

{*******************************************************}

uses Crt;

const

c:array[1..4] of real = (0,0.5,0.5,1);

type

coef = array[0..4] of real;

var

i,j,m:integer;

a,b,h,x,y,y1,y2,y3:real;

k0,k:coef;

ch:char;

{-----------------------SUBROUTINES---------------------}

{ Y = F ( x ,y ) (f = x+y) }

function f(x,y:real):real;

BEGIN

f := x + y

END;

{-------------------------------------------------------}

procedure Pausa;

BEGIN

WRITELN;WRITELN ( 'Для продолжения нажмите любую клавишу ...');

REPEAT ch := readkey UNTIL ch <> '';

END;

{------------------ОСНОВНАЯ ПРОГРАММА-------------------}

BEGIN

ClrScr;

WRITELN ( 'Введите значения концов отрезка [a,b]');

READ (a,b);

WRITELN ( 'Введите начальное значение функции y0 при x=x0 ');

READ (y);

WRITELN (' Введите число значений функции на промежутке [a,b]');

READ (m);

x:= a; h:= (b-a) / m; y1:= y; y2:= y; y3:=y;

WRITELN (' Метод Эйлера Метод Э.-Коши Метод Р.-Кутта');

WRITELN ( 'x=' ,x:5:2,' y1=',y:9:6,' y2=',y2:9:6,' y3=',y3:9:6);

FOR i:= 1 TO m DO

BEGIN

y1:= y1 + h*f(x,y1);

FOR j:=1 TO 2 DO

k0[j]:=h*f(x+2*c[j]*h , y2+2*c[j]*k0[j-1]);

y2:= y2+(k0[1]+k0[2]) / 2 ;

FOR j:=1 TO 4 DO

k[j]:= h*f(x+c[j]*h, y3+ c[j]*k[j-1]);

y3:= Y3+ (k[1]+2*k[2]+2*k[3]+k[4]) / 6;

x:= x+h;

WRITELN ('x=',x:5:2,' y1=',y1:9:6,' y2=',y2:9:6,' y3=',y3:9:6);

END;

PAUSA;

END.


Библиографический список

 

1. Демидович Б.П. Основы вычислительной математики/ Б.П.Демидович, И.А. Марон. – М.: Наука, 1970. – 664 с.

2. Демидович Б.П. Численные методы анализа/ Б.П. Демидович, И.А. Марон, Э.З. Шувалова. – М: Наука, 1967. – 368 с.

3. Копченова Н.В. Вычислительная математика в примерах и задачах/ Н.В. Копченова, И.А. Марон. – М: Наука, 1972. – 367 с.

4. Рачков А.Д. Введение в численные методы/ А.Д. Рачков. – Новосибирск НЭТИ, 1992. – 47 с.

5. Ракитин В.И. Практическое руководство по методам вычислений/ В. И. Ракитин, В.Е. Первушин.- М: Высш., 1992. – 383 с.

6. Плис А.И. MATHCAD 2000/ А.И. Плис, Н.А. Сливина. – М: Финансы и

статистика, 2000. – 655 c.

7. Дьяконов В.П. MATHCAD 2000/ В.П. Дьяконов. – СПб.: Питер, 2000. – 635 с.

 

 

Последнее изменение этой страницы: 2016-06-10

lectmania.ru. Все права принадлежат авторам данных материалов. В случае нарушения авторского права напишите нам сюда...