Лабораторная работа 7 (4 часа) Численное решение задачи Коши для одного дифференциального уравнения 1 Лабораторная работа 7 (4 часа)

• Название: Лабораторная работа 7 (4 часа) Численное решение задачи Коши для одного дифференциального уравнения 1 Лабораторная работа 7 (4 часа)
• Дата: 30.04.2019
• Размер: 0.61 Mb.
• Имя файла: ode_ru.pdf
• Тип: Лабораторная работа; #75769

Лабораторная работа № 7 (4 часа)
Численное решение задачи Коши для одного дифференциального уравнения
1
Лабораторная работа № 7 (4 часа)
Численное решение задачи Коши для одного
дифференциального уравнения
Цель работы:
получение практических навыков построения алгоритмов численного решения обыкновенных дифференциальных уравнений, оценки погрешности решения, автоматического выбора шага интегрирования, программной реализации алгоритмов на компьютере, сравнение эффективности различных методов.
Краткие теоретические сведения
Задачей Коши называют задачу решения дифференциального уравнения:
( , )
dy
f t y
dt

,
(1) где
( , )
f t y
- заданная непрерывная функция двух аргументов, с начальным условием:
0
(0)
y
y

(2)
При численном решении задачи Коши по переменной t вводят сетку 0<t
1
<t
2
<...<t
n
и ищут значения неизвестной функции в узлах сетки y
1
, y
2
,..,y
n
В явном методе Эйлера приближенное решение задачи Коши находится по рекуррентной формуле:
1
( ,
)
i
i
i
i
i
y
y
h f t y

 
,
(3) где h
i
=t
i+1
–t
i
Локальная погрешность явного метода Эйлера равна:
2
( )
2
i
i
i
y t
R
h


(4)
В неявном методе Эйлера приближенное решение задачи Коши находится в результате решения нелинейного уравнения:
1 1
1
(
,
)
i
i
i
i
i
y
y
h f t
y



 
(5)
Локальная погрешность неявного метода Эйлера равна:
2 1
(
)
2
i
i
i
y t
R
h


 
(6)
В классическом методе Рунге-Кутта четвертого порядка приближенное решение задачи Коши находится по рекуррентной формуле:
1 1
2 3
4
(
2 2
)
6
i
i
i
h
y
y
k
k
k
k

 



,
(7) где
1
( ,
)
i
i
k
f t y

;

Лабораторная работа № 7 (4 часа)
Численное решение задачи Коши для одного дифференциального уравнения
2 2
1
,
2 2
i
i
i
i
h
h
k
f t
y
k









;
3 2
,
2 2
i
i
i
i
h
h
k
f t
y
k









;


4 3
,
i
i
i
i
k
f t
h y
h k



При использовании такого метода на каждом шаге вычисляются 4 значения функции f(t,y). Метод имеет четвертый порядок точности, то есть R
i
=0(h
i
5
).
В самом общем случае метод Рунге-Кутта порядка точности p строится по рекуррентной формуле:
1 1
m
i
i
i
j
j
j
y
y
h
c k


 

,
(8) где
1 1
,
j
j
i
j i
i
i
jl l
l
k
f t
a h y
h
b k












,
(9)
a
1
=0, a
j и c
j
— константы, b
jl
— элементы нижнетреугольной матрицы, такой, что каждое
k
j
получается по предыдущим значениям k
l
Формулы (8) и (9) содержат


1 2
2 1



m
m
m
коэффициентов, подлежащих определению:


2 1

m
m
коэффициентов b
jl
, m–1 коэффициентов a
j
, и m коэффициентов c
j
Чтобы определить эти коэффициенты все функции k
j
раскладывают в ряд Тейлора в окрестности точки (y
i
, t
i
). Эти разложения подставляют в формулу (8) и результат сравнивают с рядом Тейлора для функции y(t
i+1
). Поскольку локальная погрешность определяется как разность:
R
i
=y(t
i+1
)-y
i+1
, то ставится требование, чтобы коэффициенты при всех
p
l
h
l
i
,
0
,

в разложениях для
y(t
i+1
) и y
i+1
были равны. Это требование дает систему уравнений относительно коэффициентов b
jl
, a
j
и c
j
. При этом соответствующие системы уравнений относительно
b
jl
, a
j
и c
j
для p=2,3,4 могут быть решены соответственно для m=2,3,4. Следовательно, что для того, чтобы построить метод порядка точности p для p=2,3 и 4 достаточно m=p обращений к функции f(y,t). Однако для p=5 эта система может быть решена только для
m

6. Следовательно, для метода Рунге-Кутта 5 порядка точности требуется, по крайней мере, 6 вызовов функции на каждом шаге. Аналогично для p=6 m

7.
Возможны различные способы выбора коэффициентов, приводящие к методам
Рунге-Кутта пятого порядка точности. Особый интерес представляют коэффициенты, приведенные в таблице 7.1, которые были впервые найдены Фелбергом, поскольку они позволяют при одном и том же выборе k
j
, но других с
j
* построить метод четвертого порядка точности. Как будет показано ниже, это дает возможность апостериорно оценить локальную погрешность решения. Позже были найдены и другие наборы коэффициентов,

Лабораторная работа № 7 (4 часа)
Численное решение задачи Коши для одного дифференциального уравнения
3 обладающие такой же особенностью (см. например таблицу 7.2). Многие современные алгоритмы строятся на этом же принципе, но для m=7, имеющие более высокую устойчивость.
Таблица 7.1. Коэффициенты Фелберга (Fehlberg) методов Рунге-Кутта 4 и 5 порядков точности.
j
a
i
b
jl
c
j
c
j
*
l
1 2
3 4
5 1
0 135 16 216 25 2
4 1
4 1
0 0
3 8
3 32 3
32 9
12825 6656 2565 1408 4
13 12 2197 1932 2197 7200

2197 7296 56430 28561 4104 2197 5
1 216 439
-8 513 3680 4104 845

50 9

5 1

6 2
1 27 8

2 2565 3544

4104 1859 40 11

55 2
0
Таблица 7.2. Коэффициенты Кеша-Карпа (Cash-Karp) методов Рунге-Кутта 4 и 5 порядков точности.
j
a
i
b
jl
c
j
c
j
*
l
1 2
3 4
5 1
0 378 37 27648 2825 2
5 1
5 1
0 0
3 10 3
40 3
40 9
621 250 48384 18575 4
5 3
10 3
10 9

5 6
594 125 55296 13525 5
1 54 11

2 5
27 70

27 35 0
14336 277 6
8 7
55296 1631 512 175 13824 575 110592 44275 4096 253 1771 512 4
1

Лабораторная работа № 7 (4 часа)
Численное решение задачи Коши для одного дифференциального уравнения
4
При использовании какого-либо из методов интегрирования задачи Коши задается точность, требуемая в решении. Однако, как правило, нет информации о величине шага, необходимой, чтобы достигнуть ее. Более того, требования к величине шага обычно меняются в процессе решения. Поэтому существенно, чтобы алгоритм включал автоматический выбор шага. Для этого необходимо апостериорно, т.е. после проведения вычисления, оценить погрешность полученного результата. Существует несколько подходов к практической оценке погрешности численного интегрирования дифференциальных уравнений.
Простейший подход основан на применении первого правила Рунге. Для этого на каждом шаге при переходе от t
i
к t
i+1
вычисления проводятся дважды с шагами h и
h
q
Полученные значения y
i
h

1,
, y
i
h
q

1,
служат для сравнения достигнутой точности на этом шаге:
1,
1,
1
h
i
h
i
q
i
p
y
y
R
q





,
(10) где p – порядок точности метода.
При использовании неявных методов интегрирования дифференциальных уравнений применение правила Рунге (10) приводит к большим машинным затратам, поскольку при измельчении сетки приходится в q раз больше решать нелинейных уравнений. Поэтому в этих случаях стараются применять другие апостериорные оценки погрешности. Например, при использовании неявного метода Эйлера учитывают, что если шаг h
i
мал, то можно считать, что

 

y
t
y
t
i
i
( )
(
)
1
. Следовательно, как вытекает из формул (4) и (6), погрешности явного и неявного методов почти равны по абсолютной величине, но противоположны по знаку. Поэтому после получения решения y
i+1
неявного метода Эйлера можно вычислить значение u
i+1
по явной схеме и оценить погрешность как:
R
u
y
i
i
i




(
)
1 1
2
(11)
Другой способ оценки погрешности состоит в получении оценки производной, входящей в выражение для погрешности. Так, например, для метода Эйлера необходимо оценить

y
t
i
( )
. Для этого можно воспользоваться конечно-разностной аппроксимацией производной функции f(y,t):
 
 






,
,
i
i
i
i
i
i
i
i
i
i
t t
t t
t t
t t
dy t
d
d
f dy
f
f
f
y t
f t y
f t y
dt
dt
dt
y dt
t
y
t



























, где

 

i
i
i
i
i
i
t
t
y
y
t
f
y
y
t
f
y
f
i








,
,
,

 

i
i
i
i
i
i
t
t
t
y
t
f
y
t
t
f
t
f
i








,
,
Шаг конечных разностей целесообразно выбирать так, чтобы внести возмущение примерно в половину разрядов мантиссы:
;
i
маш
i
i
маш
i
t
t
y
y







Лабораторная работа № 7 (4 часа)
Численное решение задачи Коши для одного дифференциального уравнения
5
Другой способ оценки

y
t
i
( )
состоит в использовании полученных значений y
i-1
,
y
i
, y
i+1
. Дважды дифференцируя интерполяционный полином Лагранжа, проходящий через точки (y
i-1
,t
i-1
), (y
i
, t
i
), (y
i+1
,t
i+1
) получим


















)
(
)
(
2
)
(
1 1
1 1
1 1
i
i
i
i
i
i
i
i
i
i
i
i
h
h
h
y
h
h
y
h
h
h
y
t
y
Иной способ оценки погрешности состоит в одновременном использовании двух методов разного порядка точности. Предположим, что используются два метода с порядками точности р и р+1 соответственно. Можно показать, что в этом случае локальная погрешность решения задачи Коши может быть оценена как:
1
*
1




i
i
i
y
y
R
,
(12) где
1
*
1
,


i
i
y
y
— решения, полученные методами с порядками точности р и р+1 соответственно.
Следует иметь в виду, что если используемые методы различных порядков точности используют различные сетки узлов на элементарном отрезке, то применение
(12) не имеет преимуществ по сравнение с правилом Рунге поскольку существенно не уменьшает числа обращений к функции f(y,t). Поэтому важно использовать такие методы различных порядков, которые построены на одной и той же сетке узлов. Так, например, удобно использовать методы Рунге-Кутта 4 и 5 порядка, построенным по формулам:


 
6 6
6 5
5 4
4 3
3 2
2 1
1 1
0
i
i
i
i
h
k
c
k
c
k
c
k
c
k
c
k
c
h
y
y









,
(13) где
1
( ,
)
i
i
k
f t y

;


1 21 2
2
,
k
b
h
y
h
a
t
f
k
i
i
i
i



;




2 32 1
31 3
3
,
k
b
k
b
h
y
h
a
t
f
k
i
i
i
i




;




3 43 2
42 1
41 4
4
,
k
b
k
b
k
b
h
y
h
a
t
f
k
i
i
i
i





;




4 54 3
53 2
52 1
51 5
5
,
k
b
k
b
k
b
k
b
h
y
h
a
t
f
k
i
i
i
i






;




5 65 4
64 3
63 2
62 1
61 6
6
,
k
b
k
b
k
b
k
b
k
b
h
y
h
a
t
f
k
i
i
i
i







Коэффициенты b
jl
, a
j
и c
j
для методов Фелберга и Кеша-Карпа приведены в таблицах 1 и 2.
Как было показано выше, для вышеприведенных функций
6
,
1
,

j
k
j
могут быть найдены коэффициенты
6
,
1
,
*

j
c
j
приводящие к методу четвертого порядка, т.е.:

  
5 6
*
6 5
*
5 4
*
4 3
*
3 2
*
2 1
*
1
*
1 0
i
i
i
i
h
k
c
k
c
k
c
k
c
k
c
k
c
h
y
y









Следует отметить, что на практике значение
*
1

i
y
не вычисляется, а погрешность на каждом шаге, как это вытекает из (12), оценивают из выражения:






6 1
*
j
j
j
j
i
i
k
c
c
h
R
(14)

Лабораторная работа № 7 (4 часа)
Численное решение задачи Коши для одного дифференциального уравнения
6
После оценки погрешности по одному из вышеприведенных способов, проверяется выполнение условия:
1



i
i
y
R
(15)
Если оно не выполняется, то шаг уменьшают в два раза или в
1 1



p
i
i
R
y
раз. Уменьшение шага выполняют до тех пор, пока не выполнится условие (15). После чего переходят к новому шагу. При этом новый шаг h
i+1
выбирают по правилу зон:



















,
,
;
,
1 1
1 1
1
p
i
i
i
i
i
p
i
i
q
y
R
qh
y
R
q
h
h
(16) или из соотношения:
1 1
1





p
i
i
i
i
R
y
h
h
(17)
При прогнозировании шага по правилу (17) шаг ближе к оптимальному, но общее число возвратов обычно больше, чем при использовании правила зон (16).
Следует отметить, что для интегрирования дифференциальных уравнений с быстро меняющимися функциями f(t,y) следует ограничивать шаг интегрирования заданной величиной h
max
Рабочее задание
1.
Построить алгоритмы численного решения задачи Коши явным, неявным методом
Эйлера и методом Рунге-Кутта с автоматическим выбором шага интегрирования для достижения заданной локальной относительной погрешности

либо локальной абсолютной погрешности

. Апостериорную оценку погрешности для явного метода
Эйлера провести по правилу Рунге либо по оценке производной, входящей в выражение для погрешности. Для неявного метода Эйлера погрешность следует оценивать по формуле (11), а для метода Рунге-Кутта по формуле (14).
2.
Составить рабочую программу с использованием универсальных функций численного решения задачи Коши. Минимальный набор параметров функции решения задачи
Коши должен включать: имя функции, вычисляющей производную неизвестной функции (правую часть дифференциального уравнения); интервал времени, на котором нужно найти решение; начальное условие; погрешность решения; массив временных отсчетов, в которых найдено решение; массив решений; максимальное число временных отсчетов; количество фактически полученных временных отсчетов, в которых найдены решения с заданной погрешностью.

Лабораторная работа № 7 (4 часа)
Численное решение задачи Коши для одного дифференциального уравнения
7 3.
Набрать и отладить программу на компьютере. Отладку осуществить с использованием в качестве f(t,y) функции, для которой известно аналитическое решение.
4.
Выбрать задачу в соответствии с вариантом. Для заданной схемы сформулировать задачу Коши.
5.
Провести решение задачи Коши, полученной в пункте 4 рабочего задания, каждым методом с подсчетом числа элементарных отрезков, необходимых для достижения заданной погрешности, и требуемого числа обращений к функции f(t,y).
6.
Сравнить различные методы по эффективности.
Содержание отчета
1.
Название работы
2.
Цель
3.
Рабочее задание
4.
Математические формулировки алгоритмов решения задачи Коши.
5.
Решаемое уравнение.
6.
Текст рабочей программы.
7.
Результаты расчетов с подсчетом числа необходимых разбиений заданного отрезка на элементарные и количества вызовов функции f(y,t) для каждого из методов, график решения и входного воздействия в одной системе координат.
8.
Выводы.
Задачи
1.
Найти напряжение u
вых
в течение времени


end
t
t
,
0

для схемы, приведенной на рисунке, если переключатель S был разомкнут в течение времени


0
,



t
, а в момент времени t = 0 замыкается.
Учесть, что ток i
d
, протекающий через диод связан с напряжением u
d
зависимостью:















1 0
T
d
d
m
u
i
i
exp
,

Лабораторная работа № 7 (4 часа)
Численное решение задачи Коши для одного дифференциального уравнения
8 где i
0
— обратный ток диода,

T
— тепловой потенциал, m — коэффициент неидеальности диода. Зависимость e(t) имеет вид:
)
2
sin(
)
(
0




ft
E
t
e
, где E
0
— амплитуда, f — частота,

— фаза напряжения источника.
Значения параметров напряжения источника E
0
, f и

, сопротивления резистора R, емкость конденсатора C, параметры диода и время t
end
выбрать из таблицы в соответствии с вариантом:
№ варианта
E
0
,
В
f,
10 3
Гц

R,
10 3
Ом
C,
10
-6
Ф
i
0
,
10
-12
А
m

T
,
10
-3
В
t
end
,
10
-3
c
1 5
1

/4 5
1 3000 1.7 26 3
2 10 6
6

/5 1.8 0.5 1
1.8 26 0.5 3
15 3
-

/8 0.8 2
10 1.6 26 1
4 3
2

/3 0.2 5
100 1.9 26 1.5 2.
Найти напряжение на диоде u
d
для схемы, приведенной на рисунке. Для этого заменить диод его эквивалентной схемой, приведенной на рисунке.
Учесть, что ток источника j описывается уравнением:
0
exp
1
pn
T
u
j
i
m

















, где i
0
— обратный ток диода,

T
— тепловой потенциал, m — коэффициент неидеальности диода.
Сопротивление базы диода R
b
зависит от тока:
Iv
j
R
R
b
b


1 0
, где R
b0
— немодулированное сопротивление базы диода, I
v
— ток перехода от низких к высоким уровням инжекции.
Зависимость входного воздействия e(t) показана на рисунке.

Лабораторная работа № 7 (4 часа)
Численное решение задачи Коши для одного дифференциального уравнения
9
Значения параметров напряжения источника E
0
, Tp, сопротивления резистора R, барьерная емкость диода C, параметры диода и время t
end
выбрать из таблицы в соответствии с вариантом:
№ варианта
E
0
,
В
T
p
,
10
-6
c
R,
10 3
Ом
C,
10
-9
Ф
i
0
,
10
-12
А
m

T
,
10
-3
В
R
b0
,
10 3
Ом
I
v
,
10
-3
А
t
end
,
10
-6
c
1 5
1 1
2 5000 1.7 26 1
0.3 5
2 10 2
4 1.6 300 1.8 26 3
0.1 10 3
15 1.5 6
2 100 1.6 26 0.8 0.5 9
4 8
1 1
2.5 50 1.9 26 1.8 0.1 6
3.
Найти напряжение u
вых
в течение времени


end
t
t
,
0

для схемы, приведенной на рисунке, если переключатель S был разомкнут в течение времени


0
,



t
, а в момент времени t=0 замыкается.
Учесть, что ток i
d
, протекающий через диод связан с напряжением u
d
зависимостью:
3 2
g
,
0 0,
0
d
d
d
d
u
u
i
u


 



, где g — первеанс диода. Зависимость e(t) имеет вид:
)
2
sin(
)
(
0




ft
E
t
e
, где E
0
— амплитуда, f — частота,

- фаза напряжения источника.
Значения параметров напряжения источника E
0
, f и

, сопротивления резистора R, емкость конденсатора C, параметры диода и время t
end
выбрать из таблицы в соответствии с вариантом:
№ варианта
E
0
,
В
f,
10 3
Гц

R,
10 3
Ом
C,
10
-6
Ф
g,
10
-3
A/В
3/2
t
end
,
10
-3
c
1 50 1

/4 4
1 0.8 4
2 100 5

/3 3
0.4 1
0.6 3
150 3
-

/8 0.8 2
3 1
4 200 2
-

/4 1.2 4
8 1.5

Лабораторная работа № 7 (4 часа)
Численное решение задачи Коши для одного дифференциального уравнения
10 4. Найти напряжение u
вых
в течение времени


end
t
t
,
0

для схемы, приведенной на рисунке, если переключатель S был разомкнут в течение времени


0
,



t
, а в момент времени t=0 замыкается.
Учесть, что токи эмиттера i
e
и коллектора i
с
связаны с напряжениями u
be
и u
bc
зависимостями:

















1
exp
0
T
e
be
e
e
m
u
i
i
,

















1
exp
0
T
c
bc
c
c
m
u
i
i
, где i
e0
— обратный ток эмиттерного перехода транзистора,

T
— тепловой потенциал,
m
e
— коэффициент неидеальности эмиттерного перехода транзистора, i
c0
— обратный ток коллекторного перехода транзистора, m
c
— коэффициент неидеальности коллекторного перехода транзистора. Зависимость e(t) имеет вид:
)
2
sin(
)
(
0




ft
E
t
e
, где E
0
— амплитуда, f — частота,

— фаза напряжения источника.
Значения параметров напряжения источника E
0
, f и

, сопротивления резистора R, емкость конденсатора C, параметры транзистора и время t
end
выбрать из таблицы в соответствии с вариантом:
№ варианта
E
0
,
В
f,
10 3
Гц

R,
10 3
Ом
C,
10
-6
Ф
I
e0
,
10
-12
А
m
e
I
c0
,
10
-12
А
m
c

T
,
10
-3
В
t
end
,
10
-3
c
1 5
1

/8 4
1 1
1.7 6
1.8 26 3
2 10 5
-

/6 1
0.8 5
1.6 25 1.8 26 0.6 3
15 4
-

/8 1.5 0.5 10 1.7 70 1.9 26 0.75 4
3 2

/4 0.2 5
20 1.9 80 1.8 26 1.5

Лабораторная работа № 7 (4 часа)
Численное решение задачи Коши для одного дифференциального уравнения
11
Варианты заданий
Номер бригады
Номер задачи
Номер варианта задачи
1 1
1 2
2 1
3 3
1 4
4 1
5 1
2 6
2 2
7 3
2 8
4 2
9 1
3 10 2
3 11 3
3 12 4
3 13 1
4 14 2
4 15 3
4 16 4
4