Файл: Методы численного решения дифференциальных уравнений первого порядка.doc
Добавлен: 29.10.2018
Просмотров: 822
Скачиваний: 9
СОДЕРЖАНИЕ
Методы численного решения дифференциальных уравнений первого порядка
Метод Рунге-Кутты четвертого порядка для решения уравнения первого порядка
Методы численного решения систем дифференциальных уравнений
Периодический процесс с затуханием и с внешней силой
Решение системы ДУ методом Эйлера
Решение системы ДУ методом Кранка-Николсона
Метод Рунге-Кутта четвертого порядка для решения системы уравнений первого порядка
(18)
Таким образом, дифференциальное уравнение, описывающее затухающие колебания, имеет вид:
(19)
Это дифференциальное уравнение легко представить в виде системы дифференциальных уравнений первого порядка:
(20)
Матричная форма дифференциального уравнения:
(21)
Система дифференциальных уравнений, описывающая затухающие колебания также относится к линейным однородным дифференциальным уравнениям, имеет тривиальное решение и, более того, любое решение этой системы стремится к тривиальному – то есть к состоянию равновесия.
Периодический процесс с затуханием и с внешней силой
Уравнение, описывающее периодический процесс с затуханием и с внешней возбуждающей силой.
(22)
Иногда, если внешняя сила F(t) периодическая, такую систему называют системой вынужденных колебаний.
Запишем это дифференциальное уравнение в виде системы дифференциальных уравнений первого порядка:
(23)
Запись системы дифференциальных уравнений в матричном виде:
(24)
Эта система
дифференциальных уравнений уже не
является однородной. В правой части
имеется свободное, зависящее от аргумента,
слагаемое – столбец
,
представляющий собой вектор возбуждающей
силы. Неоднородные линейные системы
дифференциальных уравнений уже не
обязаны иметь тривиальное решение.
Более того, по истечении достаточно
большого промежутка времени, когда
начальное возмущение затухнет, останется
решение, определяющееся возбуждающей
силой.
Для решения систем дифференциальных уравнений можно применять такие же методы, как и для дифференциальных уравнений первого порядка, только при этом, необходимо модифицировать расчетные формулы с учетом того, что искомая функция является вектором, и правая часть дифференциального уравнения представляет собой векторное выражение.
Решение системы ДУ методом Эйлера
Рассмотрим систему дифференциальных уравнений вида (17):
,
где Y – вектор решений, A – матрица. В общем случае коэффициенты матрицы могут зависеть от t.
Расчетная формула по методу Эйлера для такого уравнения имеет вид:
(25)
Описание алгоритма решения системы ДУ методом Эйлера:
-
Задаем размерность системы N.
-
Объявляем вектор решений Y(N).
-
Объявляем матрицу А(N,N).
-
Объявляем вспомогательный вектор P(N).
-
Задаем начальное t0 и конечное tm значения отрезка интегрирования, массив начальных значений Y0, т.е. значения системы Y(t0) ….. Y(tm) в точке t0, а также количество шагов интегрирования m
-
Вычисляем величину шага интегрирования h.
-
Цикл по t от t0 до tm c шагом h
-
Цикл по i от 1 до N
-
Вызвать процедуру Получить_Матрицу(А,t)
-
Задать P(i)=0.
-
Цикл по j от 1 до N.
-
Вычислить P(i)=P(i)+A(i,j)*Y(t).
-
Конец цикла по j.
-
Конец цикла по i.
-
Цикл по i от 1 до N.
-
Вычислить Y(i)=Y(i)+h*P(i).
-
Конец цикла по i.
-
Вызвать процедуру Печать(Y,t)
-
Конец цикла по t.
Алгоритм решения системы ДУ методом Эйлера на естественном языке:
N=размерность_системы
Объявить Массив Y(N),P(N),A(N,N)
t0=0
tm=100
m=1000
H = (tm-t0)/m
Получить_вектор (Y,t0)
Цикл по t от t0 до tm с шагом 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
КонецПроцедуры
Метод Рунге-Кутта четвертого порядка для решения системы уравнений первого порядка
Система обыкновенных дифференциальных уравнений высшего порядка путем введения новых неизвестных может быть сведена к системе уравнений первого порядка. Рассмотрим такую систему
yi ' = Fi (x, y)
где i = 1, ..., n;
y = (y1 , y2 , ..., yn );
y(x0 ) = y (0) =(y (0)1 , ..., y (0)n ).
Методе Рунге-Кутта вычисляет значение yi k по формуле:
yi(k+1) = yi (k) + (k1(i)+2*k2(i)+2*k3(i)+k4(i))/6
где:
k1(i) = Fi (x (k), y (k))*h
k2(i) = Fi (x(k)+h/2, y(k)+k1(i)/2)*h
k3(i) = Fi (x(k)+h/2, y(k)+k2(i)/2)*h
k4(i) = Fi (x(k)+h, y(k)+k3(i))*h,
h = (xn -x0)/m
n – количество уравнений системы, m - количество шагов интегрирования. Процедура использует набор функций F(i, x, y), которые соответствуют функциям Fi (x, y) описанным выше.
Описание алгоритма метода Рунге-Кутта решения систем дифференциальных уравнений первого порядка на естественном языке:
-
Задаем начальное x0 и конечное xn значения отрезка интегрирования, массив начальных значений y0, т.е. значения системы y0(i)… y0(n) в точке x0, а также количество шагов интегрирования m.
-
Вычисляем шаг интегрирования h.
-
Задаем цикл по х от x0 до xn с шагом h.
-
На основе известных x0 и y0(1)… y0(n) вычисляем правые части системы для всех n уравнений.
-
Вычисляем коэффициенты k1 и для всех n уравнений.
-
Шаги 4 и 5 проделываются для вычисления коэффициентов k2, k3, k4 для всех n уравнений системы в соответствии с формулами метода.
-
Вычисляем y1… yт для первого шага интегрирования в соответствии с главной формулой метода.
-
Шаги 4, 5, 6, 7 повторяются для вычисления состояний системы на каждом шаге интегрирования.
-
Конец цикла по х.
-
Конец алгоритма.
Алгоритм метода Рунге-Кутта решения систем дифференциальных уравнений первого порядка на естественном языке:
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