(18)

Таким образом, дифференциальное уравнение, описывающее затухающие колебания, имеет вид:

(19)

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

(20)

Матричная форма дифференциального уравнения:

(21)

Система дифференциальных уравнений, описывающая затухающие колебания также относится к линейным однородным дифференциальным уравнениям, имеет тривиальное решение и, более того, любое решение этой системы стремится к тривиальному – то есть к состоянию равновесия.

Периодический процесс с затуханием и с внешней силой

Уравнение, описывающее периодический процесс с затуханием и с внешней возбуждающей силой.

(22)

Иногда, если внешняя сила F(t) периодическая, такую систему называют системой вынужденных колебаний.

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

(23)

Запись системы дифференциальных уравнений в матричном виде:

(24)

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

Для решения систем дифференциальных уравнений можно применять такие же методы, как и для дифференциальных уравнений первого порядка, только при этом, необходимо модифицировать расчетные формулы с учетом того, что искомая функция является вектором, и правая часть дифференциального уравнения представляет собой векторное выражение.

Решение системы ДУ методом Эйлера

Рассмотрим систему дифференциальных уравнений вида (17):

,

где Y – вектор решений, A – матрица. В общем случае коэффициенты матрицы могут зависеть от t.

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

(25)

Описание алгоритма решения системы ДУ методом Эйлера:

1. Задаем размерность системы N.

2. Объявляем вектор решений Y(N).

3. Объявляем матрицу А(N,N).

4. Объявляем вспомогательный вектор P(N).

5. Задаем начальное t₀ и конечное t_m значения отрезка интегрирования, массив начальных значений Y₀, т.е. значения системы Y(t₀) ….. Y(t_m) в точке t₀, а также количество шагов интегрирования m

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

7. Цикл по t от t₀ до t_m c шагом h

8. Цикл по i от 1 до N

9. Вызвать процедуру Получить_Матрицу(А,t)

10. Задать P(i)=0.

11. Цикл по j от 1 до N.

12. Вычислить P(i)=P(i)+A(i,j)*Y(t).

13. Конец цикла по j.

14. Конец цикла по i.

15. Цикл по i от 1 до N.

16. Вычислить Y(i)=Y(i)+h*P(i).

17. Конец цикла по i.

18. Вызвать процедуру Печать(Y,t)

19. Конец цикла по t.

Алгоритм решения системы ДУ методом Эйлера на естественном языке:

N=размерность_системы

Объявить Массив Y(N),P(N),A(N,N)

t₀=0

t_m=100

m=1000

H = (t_m-t₀)/m

Получить_вектор (Y,t₀)

Цикл по t от t₀ до t_m с шагом h

Цикл по i=1 до N

Получить_матрицу (A,t)

P(i)=0

Цикл по j=1 до N

P(i) = P(i) + A(i,j)*Y(j)

КонецЦикла

КонецЦикла

Цикл по i=1 до N

Y(i) = Y(i) + h*P(i)

КонецЦикла

Печать_Вектор(Y, t)

КонецЦикла

Пример реализации алгоритма решения системы дифференциальных уравнений методом Эйлера на VFP:

SET DECIMALS TO 8

CLEAR

N=2

DIMENSION Y(N),P(N),A(N,N)

t0=0

tm=10

m=100

H =(tm-t0)/m

Получить_вектор(@Y,N,t0)

FOR t=t0 to tm STEP h

FOR i=1 TO N

Получить_матрицу(@A,N,t)

P(i)=0

FOR j=1 TO N

P(i) = P(i) + A(i,j)*Y(j)

ENDFOR

ENDFOR

FOR i=1 TO N

Y(i) = Y(i) + h*P(i)

ENDFOR

Печать_Вектор(@Y, t)

ENDFOR

PROCEDURE Получить_вектор

PARAMETERS Y,N,t0

nf=FOPEN("U:\Prog\VFP\Euler_Sys\Y.txt",0)

FOR i=1 TO N

Y(i)=VAL(FGETS(nf))

ENDFOR

FCLOSE(nf)

ENDPROC

PROCEDURE Получить_матрицу

PARAMETERS A,N,t

nf=FOPEN("U:\Prog\VFP\Euler_Sys\A.txt",0)

FOR j=1 TO N

FOR i=1 TO N

A(i,j)=VAL(FGETS(nf))

ENDFOR

ENDFOR

FCLOSE(nf)

ENDPROC

PROCEDURE Печать_Вектор

PARAMETERS Y,t

?" t= " + ALLTRIM(STR(t))

FOR i=1 TO alen(y)

? " y(" + ALLTRIM(STR(i)) + ")= " + ALLTRIM(STR(y(i)))

ENDFOR

ENDPROC

Пример реализации алгоритма решения системы дифференциальных уравнений методом Эйлера на VBA:

Sub Euler_Sys()

n = 2

ReDim Y(n): ReDim P(n): ReDim A(n, n)

t0 = 0

tm = 100

m = 1000

h = (tm - t0) / m

Call Получить_вектор(Y(), n, t0)

For t = t0 To tm Step h

For i = 1 To n

Call Получить_матрицу(A(), n, t)

P(i) = 0

For j = 1 To n

P(i) = P(i) + A(i, j) * Y(j)

Next j

Next i

For i = 1 To n

Y(i) = Y(i) + h * P(i)

Next i

Call Печать_Вектор(Y(), t)

Next t

End Sub

Sub Получить_вектор(Y(), n, t0)

Open "U:\Prog\VFP\Euler_Sys\Y.txt" For Input As #nf

For i = 1 To n

Input #nf, Y(i)

Next i

Close #nf

End Sub

Sub Получить_матрицу(A(), n, t)

Open "U:\Prog\VFP\Euler_Sys\A.txt" For Input As #nf

For j = 1 To n

For i = 1 To n

Input #nf, A(i, j)

Next i

Next j

Close #nf

End Sub

Sub Печать_Вектор(Y(), t)

Debug.Print " t= " & CStr(t)

For i = 1 To UBound(Y)

Debug.Print " y(" & CStr(i) & ")= " & CStr(Y(i))

Next i

End Sub

Решение системы ДУ методом Кранка-Николсона

Рассмотрим систему дифференциальных уравнений вида (17):

,

где Y – вектор решений, A – матрица. В общем случае коэффициенты матрицы могут зависеть от t.

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

1) ,

2) ,

3) .

В результате проделаных упрощений получаем расчетную формулу для решения системы ДУ методом Кранка-Николсона:

(26)

Алгоритм решения системы ДУ методом Кранка-Николлсона на естественном языке:

N=значение_размерности_системы_ДУ

Массив A(N,N), Y(N), E(N,N),B1(N,N),B2(N,N),B3(N,N)

Получить_единичную_матрицу(@E,N)

Задать_значения_матрице(@B1,0,N)

Задать_значения_матрице(@B2,0,N)

Задать_значения_матрице(@B3,0,N)

A=левая_граница

B=правая_граница

H=значение_шага

X=a

Получить_начальные_значения(@Y,N,a)

Цикл по x от a до b с шагом h

Получить_матрицу_правых_частей(@A,x,N)

Цикл по i=1 до N

Цикл по j=1 до N

B1(i,j)=E(i,j)-0.5*h*A(i,j)

КонецЦикла

КонецЦикла

Обратить_матрицу(@B1,N)

Цикл по i=1 до N

Цикл по j=1 до N

B2(i,j)=E(i,j)+0.5*h*A(i,j)

КонецЦикла

КонецЦикла

Умножить_матрицы(@B1,@B2,@B3,N)

Умножить_матрицу_на_вектор(@B3,@Y,@Y1,N)

Печать_решения(@Y1,x+h,N)

КонецЦикла

Процедура Получить_единичную_матрицу(E,N)

Локальные переменные i,j

Цикл по i=1 до N

Цикл по j=1 до N

Если i=j тогда

E(i,j)=1

Иначе

E(i,j)=0

КонецЕсли

КонецЦикла

КонецЦикла

КонецПроцедуры

Процедура Задать_значения_матрице(B,X,N)

Локальные переменные i,j

Цикл по i от 1 до N

Цикл по j от 1 до N

B(i,j)=X

КонецЦикла

КонецЦикла

Процедура Получить_начальные_значения(Y,N,a)

Y(1)=нач.значение_1

Y(2)=нач.значение_2

…

Y(N)=нач.значение_N

КонецПроцедуры

Процедура Получить_матрицу_правых_частей(A,x,N)

A(1,1)= значение 1_1 как функция от x

A(1,2)= значение 1_2 как функция от x

…

A(N,N)= значение N_N как функция от x

КонецПроцедуры

Процедура Обратить_матрицу(B1,N)

*** выполняется обращение матрицы B1

*** обращенная матрица подставляется вместо B1

КонецПроцедуры

Процедура Умножить_матрицы(A,B,C,N)

Локальные переменные i,j,k

Цикл по i от 1 до N

Цикл по j от 1 до N

C(i,j)=0

Цикл по k от 1 до N

C(i,j) = C(i,j)+A(i,k)*B(k,j)

КонецЦикла

КонецЦикла

КонецЦикла

КонецПроцедуры

Процедура Умножить_матрицу_на_вектор(A,B,C,N)

Локальные переменные i,j,k

Цикл по i от 1 до N

C(i)=0

Цикл по j от 1 до N

C(i)= C(i)+A(i,j)*B(j)

КонецЦикла

КонецЦикла

КонецПроцедуры

Процедура Печать_решения(Y,x,N)

*** выполняется вывод значений x и компонент вектора Y

КонецПроцедуры

Метод Рунге-Кутта четвертого порядка для решения системы уравнений первого порядка

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

y_i ' = F_i (x, y)

где i = 1, ..., n;

y = (y1 , y2 , ..., yn );

y(x0 ) = y (0) =(y (0)1 , ..., y (0)n ).

Методе Рунге-Кутта вычисляет значение y_i ^k по формуле:

y_i^(k+1) = y_i ^(k) + (k₁⁽ⁱ⁾+2*k₂⁽ⁱ⁾+2*k₃⁽ⁱ⁾+k₄⁽ⁱ⁾)/6

где:

k₁⁽ⁱ⁾ = F_i (x ^(k), y ^(k))*h

k₂⁽ⁱ⁾ = F_i (x^(k)+h/2, y^(k)+k₁⁽ⁱ⁾/2)*h

k₃⁽ⁱ⁾ = F_i (x^(k)+h/2, y^(k)+k₂⁽ⁱ⁾/2)*h

k₄⁽ⁱ⁾ = F_i (x^(k)+h, y^(k)+k₃⁽ⁱ⁾)*h,

h = (x_n -x₀)/m

n – количество уравнений системы, m - количество шагов интегрирования. Процедура использует набор функций F(i, x, y), которые соответствуют функциям F_i (x, y) описанным выше.

Описание алгоритма метода Рунге-Кутта решения систем дифференциальных уравнений первого порядка на естественном языке:

1. Задаем начальное x₀ и конечное x_n значения отрезка интегрирования, массив начальных значений y₀, т.е. значения системы y0(i)… y0(n) в точке x₀, а также количество шагов интегрирования m.

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

3. Задаем цикл по х от x₀ до x_n с шагом h.

4. На основе известных x₀ и y0(1)… y0(n) вычисляем правые части системы для всех n уравнений.

5. Вычисляем коэффициенты k₁ и для всех n уравнений.

6. Шаги 4 и 5 проделываются для вычисления коэффициентов k₂, k_3, k₄ для всех n уравнений системы в соответствии с формулами метода.

7. Вычисляем y₁… y_т для первого шага интегрирования в соответствии с главной формулой метода.

8. Шаги 4, 5, 6, 7 повторяются для вычисления состояний системы на каждом шаге интегрирования.

9. Конец цикла по х.

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

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

N=размерность_системы

Массивы Y(N),F(N),K1(N),K2(N),K3(N),K4(N)

a = значение_левой_границы

b = значение_правой_границы

h = значение_шага

Цикл по i=1 до N

F(i)=0

K1(i)=0

K2(i)=0

K3(i)=0

K4(i)=0

Y(i) = начальное_значение_i-ой_компоненты_функции

КонецЦикла

Цикл по x от a до b с шагом h

RP(x,@F,@Y,N)

Цикл по I от 1 до N

K1(i) = F(i)*h

Y1(i)= Y(i) +K1(i) /2

КонецЦикла

x1= x +h/2

RP(x1,@F,@Y1,N)

Цикл по i от 1 до N

K2(i) = F(i)*h

Y2(i)= Y(i)+K2(i) /2

КонецЦикла

RP(x1,@F,@Y2,N)

Цикл по i от 1 до N

K3(i)=F(i)*h

Y4(i)=Y(i)+K3(i)

КонецЦикла

x4=x+h

RP(x4,@F,@Y4,N)

Цикл по i от 1 до N

K4(i) = F(i)*h

Y(i) = Y(i) +(K1(i) +2*K2(i) +2*K3(i) +K4(i) )/6

КонецЦикла

Печать_решения(@Y,N,x4)

КонецЦикла

Процедура RP(x,@F,@Y,N)

Локальные переменные i

Цикл по i=1 до N

F(i) = Выражение от Y(k), x

КонецЦикла

КонецПроцедуры

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

set decimals to 10

clear

n=100

y0=0.5

m=2*n

DIMENSION y(m)

FOR i=1 TO m STEP 2

y(i)=0

y(i+1)=1

NEXT

runge_sys(0,2*PI(),@y,m,n)

PROCEDURE runge_sys

PARAMETERS xn,xe,y0,m,n

DIMENSION y0(m)

LOCAL ARRAY k1(m),k2(m),k3(m),k4(m),y1(m),y0_(m),f(m)

LOCAL i,x0,h

h=(xe-xn)/(n-1)

FOR x0 = xn TO xe STEP h

правые_части(x0,@y0,@f,m)

FOR i=1 TO m

k1(i)=h*f(i)

y0_(i) = y0(i)+k1(i)/2

NEXT

правые_части(x0+h/2,@y0_,@f,m)

FOR i=1 TO m

k2(i)=h*f(i)

y0_(i) = y0(i)+k2(i)/2

NEXT

правые_части(x0+h/2,@y0_,@f,m)

FOR i=1 TO m

k3(i)=h*f(i)

y0_(i) = y0(i)+k3(i)

NEXT

правые_части(x0+h,@y0_,@f,m)

FOR i=1 TO m

k4(i)=h*f(i)

dk=1/6*(k1(i)+2*k2(i)+2*k3(i)+k4(i))

y1(i)=y0(i)+dk

NEXT

вывод_решения(x0+h,@y1,m)

FOR i=1 TO m

y0(i) = y1(i)

NEXT

NEXT

RETURN

PROCEDURE правые_части

LPARAMETERS x,y,f,n

DIMENSION y(n),f(n)

*? "x=",x

FOR i=1 TO 200 STEP 2

f(i)=y(i+1)

f(i+1)=-y(i)

NEXT

RETURN

ENDPROC

PROCEDURE вывод_решения

PARAMETERS x,y,m

DIMENSION y(m)

? x,STR(y(1),10,7),STR(f(x),10,7),STR(ABS(y(1)-f(x)),10,7)

RETURN

ENDPROC

FUNCTION f

PARAMETERS x

RETURN SIN(x)

ENDFUNC

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

Sub runge_sys()

Const PiNumber = 3.14159265358979

Dim y()

Dim k1(), k2(), k3(), k4(), y1(), y_tmp(), f()

dec = 8

n = 100

m = 2 * n

x1 = 0: xn = 2 * PiNumber

ReDim k1(m): ReDim k2(m): ReDim k3(m): ReDim k4(m)

ReDim y1(m): ReDim y_tmp(m): ReDim f(m)

ReDim y(m)

For i = 1 To m Step 2

y(i) = 0

y(i + 1) = 1

Next i

h = (xn - x1) / (n - 1)

For x0 = x1 To xn Step h

Call правые_части(x0, y(), f, m)

For i = 1 To m

k1(i) = h * f(i)

y_tmp(i) = y(i) + k1(i) / 2

Next i

Call правые_части(x0 + h / 2, y_tmp, f, m)

For i = 1 To m

k2(i) = h * f(i)

y_tmp(i) = y(i) + k2(i) / 2

Next i

Call правые_части(x0 + h / 2, y_tmp, f, m)

For i = 1 To m

k3(i) = h * f(i)

y_tmp(i) = y(i) + k3(i)

Next i

Call правые_части(x0 + h, y_tmp, f, m)

For i = 1 To m

k4(i) = h * f(i)

dk = 1 / 6 * (k1(i) + 2 * k2(i) + 2 * k3(i) + k4(i))

y1(i) = y(i) + dk

Next i

Call вывод_решения(x0 + h, y1, m, dec)

For i = 1 To m

y(i) = y1(i)

Next i

Next x0

End Sub

Sub правые_части(x, y, f, n)

'Debug.Print "x=", x

For i = 1 To 200 Step 2

f(i) = y(i + 1)

f(i + 1) = -y(i)

Next i

End Sub

Sub вывод_решения(x, y, m, dec)

dStr = "0."

For i = 1 To dec

dStr = dStr + "0"

Next i

Debug.Print Format(x, dStr), Format(y(1), dStr), Format(f(x), dStr), Format(Abs(y(1) - f(x)), dStr)

'Return

End Sub

Function f(x)

f = Sin(x)

End Function