Лабораторная работа №9. Численное решение дифференциальных уравнений

Методы численного решения дифференциальных уравнений первого порядка

Дифференциальным уравнением первого порядка называют уравнение вида:

. (1)

При численном решении дифференциального уравнения первого порядка вида (1) с начальным условием Y (x₀) = y₀ (задача Коши) сначала выбираем фиксированное приращение аргумента h = (x_n -x₀) /n, где x₀ – начальная точка интервала интегрирования, x_n - конечная точка интервала интегрирования, n – число шагов. Далее, используя один из итерационных методов численного интегрирования дифференциальных уравнений первого порядка, находим значение Y (x_n).

Метод Эйлера

Метод Эйлера относится к наиболее простым базовым методам решения дифференциальных уравнений. Этот метод основан на представлении производной в конечно-разностной форме:

(2)

Так, для дифференциального уравнения вида (1) при достаточно малых значениях приращения аргумента x, можно, представив производную в виде отношения приращения функции к приращению аргумента, получить следующее приближенное представление в виде конечно-разностного уравнения:

(3)

Из этого представления легко выразить «новое» значение функции (значение функции в точке x+h) – расчетную формулу метода Эйлера:

(4)

Описание алгоритма численного интегрирования дифференциального уравнения первого порядка методом Эйлера:

Задаем начальное x₀ и конечное x_n значения отрезка интегрирования, начальное y₀, т.е. значение функции y в точке x₀, а также количество шагов интегрирования n.

Вычисляем величину шага интегрирования h = (x_n -x₀) /n.

Задаем цикл по k от 2 до n+1.

В цикле, на основе значения x_k-1 вычисляем x_k, т.е. x₁ = x₀ + h, x₂ = x₁ + h и т.д.

В цикле на основе значения Y_k-1 по рекуррентной формуле вычисляем y_k, т.е. Y₁ = Y₀ + h * F (x₀,Y₀), Y₂=Y₁+h*F (x ₁ ,Y ₁) и т.д. Значение F (x_k,Y_k) должно вычисляться в функции, соответствующей правой части дифференциального уравнения. Конец цикла.

Конец алгоритма.

 

Алгоритм численного интегрирования дифференциального уравнения первого порядка методом Эйлера на естественном языке:

x0=начальное_значение_х

xn=конечное_значение_х

y0=начальное_значение_Y

n=количество_шагов_интегрирования

 

h=(xn-x0)/n

 

Определить Массив X(n+1)

Определить Массив Y(n+1)

X(1)=x0

Y(1)=y0

 

Цикл по i от 2 до n+1

X(i)=X(i-1)+h

Y(i)=Y(i-1)+h*F(X(i-1),Y(i-1))

Конец Цикла

 

Печать " Решение: Значения х Значения Y(х):"

Цикл по i от 1 до n+1

Печать X(i),Y(i)

Конец Цикла

 

Определить Функцию F

Параметры x,y

f=.... *** Задается вид функции

Вернуть f

Конец Определения Функции

 

Покажем пример реализации метода Эйлера на примере уравнения

, т.е. (5)

При аналитическом интегрировании данного уравнения получаем начальное значение функции в точке x₀ =0 и вид искомой функции:

, .

Пример реализации алгоритма численного интегрирования дифференциального уравнения первого порядка методом Эйлера на VFP для уравнения (5):

SET DECIMALS TO 8

CLEAR

* Начальное и конечное значения х

x0=0

xn=1

* Начальное значение y

y0=1

* Количество итераций

n=10

 

h=(xn-x0)/n

 

DIMENSION X(n+1),Y(n+1)

X(1)=x0

Y(1)=y0

 

*Y(i+1)=Y(i)+h*F(X(i),Y(i))

FOR i=2 TO n+1

X(i)=X(i-1)+h

Y(i)=Y(i-1)+h*F(X(i-1),Y(i-1))

ENDFOR

 

? " Решение:”

? " Значения Х Известные значения Y Вычисленные значения Y:"

FOR i=1 TO n+1

?X(i),F_test(X(i)),Y(i)

ENDFOR

 

 

function F

lPARAMETERS x,y

f=y

RETURN f

ENDFUNC

 

function F_test

lPARAMETERS x

F_test=EXP(x)

RETURN F_test

ENDFUNC

 

Пример реализации алгоритма численного интегрирования дифференциального уравнения первого порядка методом Эйлера на VBA для уравнения (5):

Sub Euler_Test()

' Начальное и конечное значения х

x0 = 0

xn = 1

' Начальное значение y

y0 = 1

' Количество итераций

n = 10

 

h = (xn - x0) / n

 

ReDim x(n + 1)

ReDim y(n + 1)

x(1) = x0

y(1) = y0

 

'Y(i+1)=Y(i)+h*F(X(i),Y(i))

For i = 2 To n + 1

x(i) = x(i - 1) + h

y(i) = y(i - 1) + h * F(x(i - 1), y(i - 1))

Next i

 

Debug.Print "РЕШЕНИЕ:"

Debug.Print "Значения Х Известные значения Y Вычисленные значения Y"

For i = 1 To n + 1

Debug.Print x(i), F_test(x(i)), y(i)

Next i

 

End Sub

 

 

Function F(x, y)

F = y

End Function

 

Function F_test(x)

F_test = Exp(x)

End Function

 

Результат работы программы:

РЕШЕНИЕ:

Значения Х Известные значения Y Вычисленные значения Y

0 1 1

0,1 1,10517091807565 1,1

0,2 1,22140275816017 1,21

0,3 1,349858807576 1,331

0,4 1,49182469764127 1,4641

0,5 1,64872127070013 1,61051

0,6 1,82211880039051 1,771561

0,7 2,01375270747048 1,9487171

0,8 2,22554092849247 2,14358881

0,9 2,45960311115695 2,357947691

1 2,71828182845905 2,5937424601

Метод Кранка-Николсона

Метод Кранка-Николсона представляет собой модификацию метода Эйлера, в которой правая часть разностной формы дифференциального уравнения представляется в виде среднего значения правой части в точке x и правой части в точке x+h:

(6)

Поскольку искомое значение функции в точке x+h присутствует и в левой и в правой части уравнения – эта формула содержит Y (x+h) неявно, поэтому, для получения расчетной формулы необходимо выразить Y (x+h) явно. Выразить явно Y (x+h) возможно не для любого вида F (Y (x) ,x). Метод Кранка-Николсона можно применять, tсли правая часть дифференциального уравнения зависит только от x, т.е. F (Y (x) ,x)= F (x).

 

Рассмотрим пример дифференциального уравнения с разрешимой правой частью, зависящей от x и Y (x):

, (7)

Запишем разностный вид этого дифференциального уравнения:

После простых преобразований:

1) ,

2) ,

3) .

Получим расчетную формулу:

. (8)

В случае, когда правая часть дифференциального уравнения зависит только от х, формула метода Кранка-Николсона имеет вид:

(9)

Например, для дифференциального уравнения с правой частью, зависящей только от x:

(10)

формула Кранка-Николсона имеет вид:

. (11)

 

Описание алгоритма численного интегрирования дифференциального уравнения первого порядка с правой частью уравнения, зависящей только от х, методом Кранка-Николсона:

Задаем начальное x₀ и конечное x_n значения отрезка интегрирования, начальное y₀, т.е. значение функции y в точке x₀, а также количество шагов интегрирования n.

Вычисляем величину шага интегрирования h.

Задаем цикл по k от 2 до n+1.

В цикле на основе значения x_k-1 вычисляем x_k, т.е. x ₁ = x₀ + h, x₂ = x ₁+ h и т.д.

В цикле на основе значения Y_{k- 1} по рекуррентной формуле вычисляем y_k, т.е. Y ₁= Y₀ + h *(F (x₀,Y₀) +F (x ₁ ,Y ₁)) /2, Y₂=Y ₁ +h* (F(x ₁ ,Y ₁) +F (x₂,Y₂)) /2 и т.д. Значение F (x_k,Y_k) должно вычисляться в функции, соответствующей правой части дифференциального уравнения. Конец цикла.

Конец алгоритма.

 

Алгоритм численного интегрирования дифференциального уравнения первого порядка с правой частью уравнения, зависящей только от х, методом Кранка-Николсона на естественном языке:

x0=начальное_значение_х

xn=конечное_значение_х

y0=начальное_значение_Y

n=количество_шагов_интегрирования

 

h=(xn-x0)/n

 

Определить Массив X(n+1)

Определить Массив Y(n+1)

X(1)=x0

Y(1)=y0

 

Цикл по i от 2 до n+1

X(i)=X(i-1)+h

Y(i)=Y(i-1)+h*(F(X(i-1),Y(i-1))+ F(X(i),Y(i)))/2

Конец Цикла

 

Печать " Решение: Значения х Значения Y(х):"

Цикл по i от 1 до n+1

Печать X(i),Y(i)

Конец Цикла

 

Определить Функцию F

Параметры x,y

f=.... *** Задается вид функции

Вернуть f

Конец Определения Функции

 

Покажем пример реализации метода Кранка-Николсона на примере уравнения (10). При аналитическом интегрировании уравнения (10) получаем начальное значение функции в точке x₀ =0 и вид искомой функции:

, . (12)

 

Пример реализации алгоритма численного интегрирования дифференциального уравнения первого порядка с правой частью, зависящей только от х, методом Кранка-Николсона на VFP для уравнения (10):

SET DECIMALS TO 8

CLEAR

* Начальное и конечное значения х

x0=0

xn=1

* Начальное значение y

y0=3

* Количество итераций

n=10

 

h=(xn-x0)/n

 

DIMENSION X(n+1),Y(n+1)

X(1)=x0

Y(1)=y0

 

*Y(i+1)=Y(i)+h*F(X(i),Y(i))

FOR i=2 TO n+1

X(i)=X(i-1)+h

Y(i)=Y(i-1)+h*F(X(i-1),Y(i-1))

ENDFOR

 

? " Решение:"

? " Значения Х Известные значения Y Вычисленные значения Y:"

FOR i=1 TO n+1

?X(i),F_test(X(i)),Y(i)

ENDFOR

 

 

function F

lPARAMETERS x,y

f=x*x+3*x

RETURN f

ENDFUNC

 

function F_test

lPARAMETERS x

F_test=2*x+3

RETURN F_test

ENDFUNC

 

Пример реализации алгоритма численного интегрирования дифференциального уравнения первого порядка с правой частью, зависящей только от х, методом Кранка-Николсона на VBA для уравнения (10):

Sub Krank_Nikolson_Test()

' Начальное и конечное значения х

x0 = 0

xn = 1

' Начальное значение y

y0 = 3

' Количество итераций

n = 10

 

h = (xn - x0) / n

 

ReDim x(n + 1)

ReDim y(n + 1)

x(1) = x0

y(1) = y0

 

'Y(x+h) = Y(x) + h*(x^2 + 3*x + (x+h)^2 + 3*(x+h))/2

For i = 2 To n + 1

x(i) = x(i - 1) + h

y(i) = y(i - 1) + h * (F(x(i - 1), y(i - 1)) + F(x(i), y(i))) / 2

Next i

 

Debug.Print "РЕШЕНИЕ:"

Debug.Print "Значения Х Известные значения Y Вычисленные значения Y"

For i = 1 To n + 1

Debug.Print x(i), F_test(x(i)), y(i)

Next i

 

End Sub

 

 

Function F(x, y)

F = x ^ 2 + 3 * x

End Function

 

Function F_test(x)

F_test = 2 * x + 3

End Function

 

Результат работы программы на VBA:

РЕШЕНИЕ:

Значения Х Известные значения Y Вычисленные значения Y

0 3 3

0,1 3,2 3,0155

0,2 3,4 3,063

0,3 3,6 3,1445

0,4 3,8 3,262

0,5 4 3,4175

0,6 4,2 3,613

0,7 4,4 3,8505

0,8 4,6 4,132

0,9 4,8 4,4595

1 5 4,835