close

Вход

Забыли?

вход по аккаунту

?

numerical systems of SSU

код для вставкиСкачать
1
2
МЕТОДЫ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ УРАВНЕНИЙ
Моделирование линейных электрических цепей основано на законах
Кирхгофа (первый закон – следствие принципа непрерывности
электрического тока; второй – следствие закона сохранения энергии) и Ома.
Целью расчетов является определения токов во всех ветвях,
напряжений в узлах, мощностей электроприемников.
Математическая модель электрической цепи представляется
эквивалентной схемой, которая содержит m ветвей и n узлов, а также ЭДС
всех источников, внутренние и внешние сопротивления.
Последовательность расчетов линейной электрической цепи с
использованием законов Кирхгофа такая:
а) для каждой ветви вводят обозначение протекающего в ней тока и
стрелками на схеме указывают условные положительные направления этих
токов;
б) для n−1 узлов составляют уравнения на основании первого закона
Кирхгофа; для одного из узлов такое уравнение не составляют, поскольку
оно является следствием уже написанных уравнений;
в) берут взаимно независимые контуры (это означает, что в каждом
новом контуре хотя бы в одной из ветвей ток не входит в ранее
рассмотренные контуры), в каждом из этих контуров выбирают условное
положительное направление обхода и обозначают его на схеме;
г) для выбранных контуров составляют уравнения по второму закону
Кирхгофа с учетом направления обхода; при правильном выборе контуров их
число должно равняться m - n + 1; при этом общее количество уравнений
должно составлять m, то есть равное числу неизвестных величин (токов);
д) решают полученную систему с m уравнений одним из методов
вычислительной математики.
Если рассчитанное значение тока в данной ветви является
положительным, то это означает, что действительное направление тока
совпадает с выбранным раньше и наоборот.
3
Согласно законам Кирхгофа, алгоритм расчетов, схемы электрической
цепи и направления обхода контуров система линейных алгебраических
уравнений (СЛАУ) имеет вид, где неизвестные − это токи в ветвях,
коэффициенты при неизвестных – параметры элементов, а свободные члены
– ЭДС источников:
В дальнейшем полученная СЛАУ решается разными методами.
3.1 Формы записи систем линейных алгебраических уравнений и
способы решения матричных уравнений
В общем виде система n линейных алгебраических уравнений с n
неизвестными записывается так:
a11x1  a12 x2  ...  a1n xn  b1,
a x  a x  ...  a x  b ,
 21 1
22 2
2n n
2


an1 x1  an2  ...  ann xn  bn .
Коротко систему (3.2) можно записать в таком виде:
(3.2)
4
n
a x
j
ij
j
 bi (i = 1, 2, …, n).
(3.3)
С использованием матричных обозначений систему (3.2) можно
записать в следующем виде:
A X  B ,
(3.4)
или в развернутом виде:
a11
a
 21


an1
a12
a22
an 2
a1n   x1  b1 
a2n   x 2 b2 
 ·   =  .
    
    
ann  xn  bn 
(3.5)
Вообще, рассматриваются три вида матричных уравнений и способы
их решения.
Матричное уравнение вида A X  B .
Для решения такого уравнения надо умножить слева обе его части на
обратную матрицу A1 : A1  A X  A1  B . При этом произведение A1  A  E
(единичная матрица), итак, E  X  A1  B , откуда
X  A1  B .
(3.6)
Матричное уравнение вида X  A  B . Надо умножить справа обе его
части на A1 : X  A A1  B  A1 . Выходит: X  E  B  A1 , откуда
X  B  A1 .
(3.7)
Матричное уравнение вида A X  B  C . При умножении обеих частей
B1
уравнения
–
слева
на
A1 ,
а
справа
на
значит
A1  A X  B  B1  A1  C  B1 , или E  X  E  A1  C  B1 , откуда
X  A1  C  B1 .
(3.8)
5
3.2 Методы решения СЛАУ
Методы разделяют на прямые, которые используют для вычисления
неизвестных конечные соотношения (формулы) и численные, которые, в
общем случае, делятся на итерационные методы и методы минимизации.
Прямые методы дают решение после выполнения заранее известного
количества операций. Эти методы простые и универсальные. Вместе с тем им
присущи следующие недостатки:
– они нуждаются в хранении в оперативной памяти сразу всей
матрицы;
– не учитывают структуру матрицы – нулевые элементы также
хранятся в памяти и над ними проводятся арифметические действия;
– происходит накопление погрешностей в процессе решения (хоть
прямые методы и называют точными), потому что вычисление на любом
этапе используют результаты предыдущих операций.
В связи с этим прямые методы используют для сравнительно
небольших (n < 200) систем с плотно заполненной матрицей и не близким к
нулю определителем.
Численные методы, как итерационные, так и минимизации – это
методы последовательных приближений. Они нуждаются в некотором
приближенном решении – начальном приближении. Дальше по
определенному алгоритму выполняется ряд итераций до получения
результата с необходимой точностью. Методы имеют преимущество перед
прямыми в следующем:
– нуждаются в хранении не всей матрицы системы, а лишь нескольких
векторов с n компонент; иногда элементы матрицы можно совсем не
сохранять, а вычислять их при необходимости;
– погрешности конечных результатов не накапливаются, потому что в
каждой итерации используется результат только предыдущей и не
используются раньше выполненные вычисления.
Но при этих преимуществах сходимость итерационных методов может
быть очень медленной.
3.2.1 Прямые метода решения СЛАУ
Прямые методы рассматриваются в вузовском курсе «Высшая
математика» и в данном курсе о них только упоминается.
6
Итак, первый из прямых методов – это метод, который использует
формулы Крамера, где неизвестные определяются как отношения
определителей.
Наиболее распространенным прямым методом решения СЛАУ на ЭВМ
является метод последовательного исключения неизвестных (метод Гаусса)
и его модификации.
Одна из разновидностей метода Гаусса (компактная схема Гаусса) –
метод LU–преобразование.
Компактные схемы Гауса используются для СЛАУ с разреженными
матрицами произвольной структуры. В этих схемах применяется
триангуляция матриц (треугольная декомпозиция) – суть которой состоит в
преобразовании матрично-векторного уравнения (1.5) у уравнения с
треугольными матрицами (методику преобразования квадратной матрицы на
произведение двух треугольных приведено ниже).
Порядок действий в методе LU–преобразование следующий:
а) выполняется замена матрицы А на произведение двух треугольных
матриц L и U; (в матрице U диагональные элементы равняются единицам);
б) решается относительно Y уравнение L Y  B ;
в) решается относительно X уравнение U  X  Y .
Методика разложения матрицы на произведение двух треугольных
матриц.
Квадратную матрицу можно разложить на две треугольных
(квадратные, в которых элементы расположенные выше (ниже) главной
диагонали, равняются нулю) и это разложение будет единственным, если
диагональным элементам одной из треугольных матриц заранее дать не
равные нулю значения (например, равные единицам).
Пусть A  T1 T2 , где Т1 и Т2 – треугольные матрицы в одной из которых
диагональными элементами являются единицы. Другие элементы матриц Т1 и
Т2 находят следующим образом:
а) перемножают матрицы Т1 и Т2 (их элементы – буквенные
обозначения);
б) приравнивают соответствующие элементы матрицы-произведения
Т1Т2 элементам матрицы А – получают уравнение: одночленные, двучленные
и т.д.;
в) решают полученные уравнения – сначала одночленные, дальше
двучленные и т.д.
7
3.2.2 Итерационные метода решения СЛАУ
Рассматриваются три итерационных метода: метод простой итерации
(последовательных приближений), метод Зейделя и метод скорейшего
спуска.
Метод простой итерации.
В общем случае задана СЛАУ (3.2), которая записана в развернутом
матричном виде (3.5). Если предположить, что диагональные элементы
матрицы А aii  0 (и = 1, 2, …, n), то можно перевести систему к
каноническому виду и потом выразить х1 через первое уравнение системы, х2
– через второе уравнение и т.д. В результате получим систему, которая
эквивалентная системе (3.2):
b a
a
a

x  1  12 x  13 x   1n x ,
 1 a
a 2 a 3
a n
11
11
11
11
b2 a21
a23
a2n

x 

x 
x  
x ,
1
3
 2 a
a22
a22
a22 n

22


a
b
a
a

x  n  n1 x  n 2 x   n,n1 x .
 n a
ann 1 ann 2
ann n1

nn
Обозначим
(3.9)
aij
bi
 βi ;  aii  αij , где i, j = 1, 2, …, n. Тогда система (3.9)
aii
запишется так:
x1  β1  α12 x2  α13 x3   α1n xn ,
x  β  α x  α x   α x ,
 2
2
21 1
23 2
2n n


 xn  β n  α n1 x1  α n2 x2   α n,n1 xn1.
(3.10)
Систему (3.10) называют приведенной к нормальному виду. Эта
система в матричной форме записи:
X  β  αX ,
или
8
α12
 x1  β1  0
 x  β   α 0
 2    2    21
    
    
 xn  βn  α nn α n 2
α1n   x1
α 2n  x2 

.
  
  
0   xn 
(3.11)
После нормализации системы проверяется условие сходимости
итерационного процесса. Признаком сходимости является условие того, что
любая из норм матрицы  меньше единицы, то есть
q  Norα1,
где q – норма матрицы α , которая может быть определена по одной из
формул:
n
q  max
 αij ,
i
j1
n
q  max  αij .
j
i1
Алгоритм метода простой итерации следующий:
– за нулевое приближение принимается столбец свободных членов:
 x1(0)  β1 
 (0)   
β
x 2    2  – нулевое приближение,
   
 (0)  β 
x   n 
 n 

дальше строятся матрицы-столбцы следующих приближений:
α
 x1(1)  β1  0
12
 (1)    α 0
2
x2   β    21
    
  β  α
α
 x(1)   n   nn n 2
 n 
α1n   x1(0) 
  
α 2n  2 x(0) – первое приближение;

  

0  x(0) 
 n 
9
α
 x1(2)  β1  0
12
 (2)    α 0
x2   β 2    21
    
  β  α
α
 x(2)   n   nn n 2
 n 

α1n   x1(1) 
  
α 2n  x2(1)  – второе приближение

  

0  x(1) 
 n 
и т.д.
Вообще, любое (k+1)-е приближения вычисляют по формуле
X (k 1)  β  αX (k ) (k = 0, 1, …, n).
(3.12)
Итерационный процесс длится до тех пор, пока не будет выполнено
условие окончания
q
max x(k 1)  x(k )  ε,
10
1 q
i
i
i
где ε – заданная абсолютная погрешность.
В методе итераций замена значений всех переменных проводится
одновременно (одновременное смещение).
Метод Зейделя.
В методе Зейделя уточненное значение х1 сразу же используется для
вычисления х2, дальше новые значения х1 и х2 используются для вычисления
х3 и т.д.
Это небольшое усовершенствование итерационной процедуры
позволяет существенно увеличить скорость сходимости.
Любое (k+1)-е приближение в методе Зейделя находится по
следующим формулам:
n
(k 1)
1
x
  1  xij(jk ) ,
j1
n
(k 1)
2
x
 2   x
(k 1)
21 1
  2 j x(kj ) ,
j2
(3.17)
n1
xn(k 1)  n  nj x(kj 1)   nn xnk ) ,
j1
где k = 0, 1, 2, …, n.
Итерации заканчиваются, когда с заданной точностью получены
одинаковые значения неизвестных в двух итерациях подряд.
Условие сходимости итерационного процесса подобно условию для
простой итерации. Итерационный процесс и его сходимость зависят от
величины элементов матрицы  следующим образом: если самая большая
сумма модулей элементов строк или самая большая сумма модулей
элементов столбцов меньше единицы, то процесс итерации для данной
системы сходится к единственному решению независимо от выбора
начального приближения.
Итак, условия сходимости можно записать так:
n
max  ij  1 (i = 1, 2, …, n) или
i
j1
n
max  ij  1 (j = 1, 2, …, n).
j
i1
11
Как и в методе простой итерации требуется привести СЛАУ к виду,
пригодному для итераций. Для выполнения условий сходимости
итерационного процесса достаточно, чтобы значение элементов αij матрицы
α при i  j были небольшие по абсолютной величине.
Это равносильно тому, что если для СЛАУ A X  B модули
диагональных коэффициентов каждого уравнения системы больше суммы
модулей всех других коэффициентов (без учета свободных членов), то
итерационные процессы для этой системы сходятся.
На примере показывается, как выполняется эквивалентное
преобразование исходной СЛАУ и получается нормализованная система в
общем случае.
Исходная СЛАУ:
4,5x1 1,8x2  3,6x3  1,7;

1,8x1  2,5x2  4,6x3  2, 2;
3,1x  2,3x 1, 2x  3,6.

1
2
3
(P)
(Q)
(R)
Выполняются следующие действия:
а) в заданной системе выделяются уравнения с коэффициентами,
модули которых большие чем сумма модулей других коэффициентов
уравнения. Каждое выделенное уравнение записывается в ту сроку новой
СЛАУ, чтобы самый большой по модулю коэффициент был диагональным. В
уравнении (Q) выполняется такое неравенство: 4,6  1,8  2,5 . Уравнение
(Q) принимается за третье уравнение новой системы.
б) остальные уравнения новой эквивалентной системы получаются
путем составления линейных независимых между собой комбинаций. Так, за
первое уравнение можно принять такую линейную комбинацию (P)+(R),
тогда:
7,6x1  0,5x2  2,4x3 1,9 .
За второе уравнение
(2Q)+(R)−(P), то есть
новой
системы
2,2x1  9,1x2  4,4x3  9,7 .
–
такую
комбинацию:
12
В результате получена преобразованная СЛАУ которая эквивалентна
исходной и удовлетворяет условиям сходимости итерационного процесса:
7,6x
1  0,5x2  2, 4x3  1,9;

2, 2x  9,1x  4, 4x  9,7;

1
2
3
1,8x  2,5x  4,6x  2, 2.

1
2
3
Для проверки этого утверждения эквивалентная система приводится к
нормальному виду и проверяется, удовлетворяется ли хоть одно из условий
сходимости.
При этом используется следующий способ: записываются
коэффициенты при неизвестных x1, x2, x3 в соответствующих уравнениях
системы как m·x, где m – число, близкое к коэффициенту при
соответствующем неизвестном и на которое легко разделить коэффициенты
при неизвестных и свободные члены (обычно 10; 100).
Например, принимается m = 10. Тогда система, приводимая к
нормальному виду, перепишется так:
10x1  2, 4x1  0,5x2  2, 4x3 1,9;

10x2  2, 2x1  0,9x2  4, 4x3  9,7;
10x  1,8x  2,5x  5, 4x  2, 2.

3
1
2
3
или
x
 1  0, 24x1  0,05x2  0, 24x3  0,19;
x  0, 22x  0,09x  0, 44x  0,97;
 2
1
2
3
x  0,18x  0, 25x  0,54x  0, 22.

3
1
2
3
Матрица α и вектор β принимают вид
0, 24 0,05 0, 24
0, 22 0,09 0, 44 ,
α  

 0,18 0, 25 0,54
0,19 
β  0,97 .


0, 22

Суммы модулей элементов строк матрицы α :
0,24 + 0,05 + 0,24 = 0,53;
0,22 + 0,09 + 0,44 = 0,75;
0,18 + 0,25 + 0,54 = 0,97.
Большая из сумм 0,97 < 1, и одно из условий сходимости
итерационного процесса выполняется. И хоть второе условие не выполняется
13
(большая сумма модулей элементов столбцов 0,24 + 0,44 + 0,54 = 1,22 > 1)
процесс итерации для рассматриваемой системы сходится к единственному
решению.
15
14
РЕШЕНИЕ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ
УРАВНЕНИЙ
7.1 Вводные замечания
Инженеру-исследователю в своей деятельности постоянно приходится
сталкиваться с дифференциальными уравнениями, так как большая часть
законов физики при их математическом моделировании сводится к
дифференциальным уравнениям. В связи с этим решение дифуравнений
является одной из важнейших математических задач.
Дифференциальным уравнением называется уравнение, в котором
искомой является функция, а уравнение задает связь между значениями
независимых переменных, искомой функцией и ее производными.
В зависимости от числа независимых переменных ДУ делятся на две
различные группы: обыкновенные ДУ, содержащие одну независимую
переменную, и уравнения в частных производных, содержащие несколько
независимых переменных.
Рассмотрим обыкновенные ДУ.
Общий вид ОДУ
F (x, y, y’, …, y(n)) = 0,
(7.1)
где x – независимая переменная;
y – функция;
y’,…, y(n) – производные;
n – наивысший порядок производной, порядок ДУ.
В ряде случаев из общей записи ДУ (7.1) удается выразить старшую
производную в явном виде, например,
y” = f (x, y, y’ ).
Такая форма записи называется уравнением, разрешенным относительно
старшей производной.
Решением ДУ (7.1) называется всякая функция y = φ(x), которая после ее
подстановки в уравнение превращает его в тождество.
Дифференциальное уравнение, как правило, имеет бесконечное
множество решений, зависящее от некоторого количества (определяемого
порядком ДУ) произвольных постоянных С1, С2, …, Сn.
Общее решение ОДУ (7.1) имеет вид
15
у = φ(x, С1, С2, …, Сn ).
Частное решение ОДУ получается из общего, если произвольным
постоянным придать определенные значения, например, для уравнения
первого порядка:
общее решение у = φ(x, С);
частное решение, если С = С0 (определенное значение) у = φ(x, С0).
Исследователя обычно интересует частное решение. Чтобы его
получить, необходимо учитывать не только особенности изменения явления,
описываемого ДУ (7.1), но и дополнительные условия, которые
характеризуют это явление. Количество дополнительных условий должно
быть равно количеству произвольных постоянных, то есть порядку n.
Дифференциальное уравнение вместе с дополнительными условиями
называется задачей. Если х является временем (что характерно для
эволюционных или динамических процессов), то необходимо учесть
состояние в начальный момент времени х0, а если х – координата, которая
может изменяться от а до b (например, прогиб балки под действием
постоянной силы), то необходимо учесть состояние системы (балки) на ее
границах – в точках а и b.
Таким образом, возможно существование двух типов дополнительных
условий и соответственно двух типов задач:
1.
Начальные условия (всего n условий) у(х0) = у0, y' (x )  y0' , …, 0
yn1(x 0)  yn1
0 .
Задача решения ОДУ (7.1) с начальными условиями,
заданными в одной точке, называется задачей Коши.
2.
Граничные (краевые) условия (всего n условий) у(а) = уа,
y' (a)  y'a, …, у(b) = уb, y' (b)  y' ,b…. Задача решения ОДУ (7.1) с граничными
условиями называется краевой или двухточечной задачей.
7.2 Как решать ДУ
Целью решения ДУ является получение функции у = φ(x). Эта функция
может быть задана графическим, аналитическим и табличным способом.
Поэтому методы решения ДУ (7.1) делятся на графические, аналитические и
численные.
Графические методы используются при качественном исследовании
ДУ и состоят в построении интегральной кривой.
16
Аналитические методы дают решение уравнения в виде формулы и
делятся на точные и приближенные.
Точные методы рассматриваются в вузовском курсе. Они дают точные
решения, но могут быть использованы для узкого класса задач.
Приближенные аналитические методы дают функцию, являющуюся
аппроксимацией точного решения. Для получения такой функции либо
заменяют исходное ДУ близким к нему, но допускающим аналитическое
решение, либо задаются видом решения.
Основной проблемой, возникающей при этом является проблема
точности. Выход состоит в использовании идеи о разбиении области
изменения независимой переменной х на малые отрезки и решении ДУ на
каждом отрезке.
Этот подход иллюстрируем следующим рисунком.
у
1
2
3
2
3
1
х
а
b
Рисунок 7.1 – Решение ДУ: 1 – точное решение уравнения; 2 –
аппроксимация решения линейной функцией для всей области; 3 –
аппроксимация решения при разбиении области изменения х[a, b] на мелкие
отрезки.
Видно, что замена линейной функции 2 кусочно-линейной функцией 3
значительно повышает точность решения.
Дальнейшее развитие рассматриваемой идеи привело к появлению
численных методов решения ДУ, результатом которых является таблица
значений функции в определенных точках (узлах).
17
7.3 Одношаговые методы решения задачи Коши
Задачу Коши можно сформулировать для ДУ первого порядка
следующим образом.
Пусть дано ДУ y’ = f (x, y) c начальным условием y(x0) = y0. Требуется
найти функцию y = f (x), удовлетворяющую как указанному уравнению, так и
начальному условию.
Обычно численное решение этой задачи получают, вычисляя сначала
значение производной в точке х0, а затем задавая малое приращение х и
переходя к новой точке x1 = x0 + h. Положение новой точки определяется по
наклону кривой, вычисленному с помощью ДУ.
Таким образом, график численного решения представляет собой
последовательности коротких прямолинейных отрезков, которыми
аппроксимируется истинная кривая y = f (x).
Сам численный метод определяет порядок действия при переходе от
данной точки кривой к следующей.
В одношаговых методах для нахождения следующей точки на кривой
y=f(x) требуется информация лишь об одном предыдущем шаге. К
одношаговым относятся метод Эйлера и методы Рунге-Кутта.
7.3.1 Метод Эйлера (схема ломаных)
Является сравнительно грубым и применяется в основном для
ориентировочных расчетов. Основан на разложении функции у в ряд Тейлора
в окрестности точки х0:
y (x0 + h) = y(x0) + h y’(x0) + 1/2 h2 y”(x0) + …
Если шаг интегрирования h мал, то члены ряда, содержащие h во второй
и более высоких степенях, малы. Поэтому ими можно пренебречь. Тогда
y (x0 + h) = y(x0) + h y’(x0).
Значение y’(x0)
yi +1 = yi + h f (xi, yi), i = 0, 1, …, n – 1.
Графически метод представляется следующим образом (рис.7.2).
18
у
уi+1
уi
х
h
xi
xi+1
Рисунок 7.2 – Графическая интерпретация метода Эйлера
Погрешность этого метода пропорциональна шагу h.
Большая погрешность – главный недостаток этого метода. Точность
метода Эйлера можно повысить, улучшив аппроксимацию производной. Это
можно сделать, используя среднее значение производной в начале и конце
интервала.
Пример 7.1.
Решить ДУ
y  2  x  y 2  x 4 на интервале [0, 1] с шагом h = 0,25
методом Эйлера. Начальные условия: y(x0) = у0 = 0.
Решение выполняется по формуле yi +1 = yi + h f (xi, yi).
Первый шаг: i = 0, x0 = 0, y0 = 0
у0 +1=1 = y0 + h (2x0 – y20 + х40) = 0 + 0,25·(2∙0 – 02 + 04) = 0.
Второй шаг: i = 1, x1 = x0 + h = 0 + 0,25 = 0,25, y1 = 0
у1 +1=2 = y1 + h (2x1 – y21 + х41) = 0 + 0,25·(2∙0,25 – 02 + 0,254) = 0,126.
Третий шаг: i = 2, x2 = x1 + h = 0,25 + 0,25 = 0,5, y2 = 0,126
19
у2 +1=3 = y2 + h (2x2 – y22 + х42) = 0,126 + 0,25·(2∙0,5 – 0,1262 + 0,54) =
0,3877.
Четвертый шаг: i = 3, x3 = x2 + h = 0,5 + 0,25 = 0,75, y3 = 0,3877
у3 +1=4 = y3 + h (2x3 – y23 + х43) = 0,3877 + 0,25·(2∙0,75 – 0,38772 + 0,754) =
0,8042.
Метод Эйлера может быть применен к решению систем
дифференциальных уравнений и ДУ высших порядков. Однако в последнем
случае ДУ должны быть приведены к системе ДУ первого порядка.
Пусть задана система двух уравнений первого порядка
 y '  f1 (x, y, z),
 '
 z  f2 (x, y, z)
с начальными условиями у(х0) = у0, z(х0) = z0.
Приближенные значения у(хi) ≈ уi и z(хi) ≈ zi находятся по формулам
yi+1 = yi + Δ yi,
zi+1 = zi + Δ zi,
где
Δ yi = hf1(xi, yi, zi), Δ zi = hf2(xi, yi, zi)
i = 0, 1, 2, …
Пример 7.2.
 y '  (z  y)x,
Применяя метод Эйлера, решить систему ДУ  '
при
z

(z

y)x

начальных условиях у(0) = 1,0000 z(0) = 1,0000 на отрезке [0, 0,6] с шагом h =
0,1. Вычисления вести с одним запасным знаком.
Покажем часть вычислений:
 первый шаг: i = 0, х0 = 0, y0 = 1,0000, z0 = 1,0000,
y0’= (z0 – y0)·x0 = (1 – 1)·0 = 0;
Δ y0 = h y0’ = 0,1·0 = 0;
z0’= (z0 + y0)·x0 = (1 + 1)·0 = 0;
Δ z0 = h z0’ = 0,1·0 = 0;
20
 второй шаг: i = 1, х1 = х0 + h = 0 + 0,1 = 0,1;
y1 = y0 + Δ y0· = 1,0000 + ·0 = 1,0000;
z1 = z0 + Δ z0· = 1,0000 + 0 = 1,0000;
y1’= (z1 – y1)·x1 = (1 – 1)·0,1 = 0;
z1’= (z1 + y1)·x1 = (1 + 1)·0,1 = 0,2;
Δ y1 = h y1’ = 0,1·0 = 0;
Δ z1 = h z1’ = 0,1·0,2 = 0,02;
 третий шаг: i = 2, х2 = х1 + h = 0,1 + 0,1 = 0,2;
y2 = y1 + Δ y1· = 1,0000 + 0 = 1,0000;
z2 = z1 + Δ z1· = 1,0000 + 0,02 = 1,02;
y2’= (z2 – y2)·x2 = (1,02 – 1)·0,2 = 0,004;
Δ y2 = h y2’ = 0,1·0,004 = 0,0004;
z2’= (z2 + y2)·x2 = (1,02 + 1)·0,2 = 0,404; Δ z2 = h z2’ = 0,1·0,404 = 0,0404.
Решение на остальных шагах приведено в табл.7.1
Таблица 7.1 – Решение системы ДУ
yi+1 =
i
xi
yi’=
yi-1+Δ yi-1 (zi – yi)·xi
0
0,1
1,0000
0
1
0,2
1,0000
0
2
0,3
1,0000
0,004
3
0,4
1,0004
0,018
4
0,5
1,0022
0,048
5
0,6
1,0070
0,1001
6
0,7
1,0170
Δ yi =
yi’·h
0
0
0,0004
0,0018
0,0048
0,0100
zi+1 =
zi’=
zi-1+Δ zi-1 (zi + yi)·xi
1,0000
0
1,0000
0,2000
1,0200
0,4040
1,0604
0,6182
1,1222
0,8498
1,2072
1,01071
1,3179
Δ zi =
zi’·h
0
0,0200
0,0404
0,0618
0,0850
0,1107
21
Пример 7.3.
Применяя метод Эйлера, составить на отрезке [1; 1,5] таблицу значений
y'
'' 
 y  0 при начальных условиях у(1) = 0,77, у’(1) = − 0,44 с шагом
ДУ y
x
h = 0,1.
С помощью подстановки у’ = z, у’’ = z’ заменим данное уравнение
системой уравнений
 y'  z,

'
z
z    y.


x
при начальных условиях у(1) = 0,77, z(1) = − 0,44.
Результаты расчета приведены в табл. 7.2
Таблица 7.1 – Решение системы ДУ
i
xi
yi
y i ’= z
Δ yi =
yi’·h
0
1
2
3
4
5
1,0
1,1
1,2
1,3
1,4
1,5
0,77
0,726
0,679
0,629
0,576
0,521
-0,44
-0,473
-0,503
-0,529
-0,551
zi
-0,044
-0,0473
-0,0503
-0,0529
-0,0551
-0,44
-0,473
-0,503
-0,529
-0,551
zi’=
- zi / xi yi
-0,33
-0,296
-0,26
-0,222
-0,183
Δ zi =
zi’·h
-0,033
-0,03
-0,026
-0,022
-0,018
7.3.2 Метод Эйлера Коши
Сначала вычисляется значение функции в следующей точке по методу
Эйлера
y*  y  hf (x , y ),
i1
i
i
i
которое затем используется для вычисления приближенного значения
*
производной в конце интервала f (xi1 , yi1
).
Вычислив среднее арифметическое между значениями производной
(наклонами ломаных) в начале и конце интервала, найдем более точное
значение уi+1 (формула Эйлера-Коши):
22
y
1
 y  h[ f (x , y )  f (x , y* )],
i1
i
i
i
i1
i1
2

где хi+1 = хi + h.
Принцип, на котором основано улучшение метода Эйлера, можно
пояснить и иначе – на основе использования разложения функции в ряд
Тейлора. Точность повышается за счет того, что сохраняются в разложении
члены, содержащие производные более высоких порядков, чем первый. Для
улучшения метода Эйлера надо знать вторую производную. Ее можно
аппроксимировать конечной разностью
y' y' (x h)  y' (x )
y (x ) 
0 0
x
h
''
0
Если ее подставить в ряд Тейлора, то результатом будет формула ЭйлераКоши.
Пример 7.4.
Решить ДУ
y  2  x  y 2  x 4 на интервале [0, 1] с шагом h = 0,25
методом Эйлера - Коши. Начальные условия: y(x0) = у0 = 0.
Решение выполняется по формулам:
1
y  y  h[ f (x , y )  f (x , y* )],
у* = y + h f (x , y ),
i +1
i
i i
i1
i
i
i
i1
i1
2


Первый шаг: i = 0, x0 = 0, y0 = 0
у 1* = y0 + h (2x0 – y20 + х40) = 0 + 0,25·(2∙0 – 02 + 04) = 0;
x1 = x0 + h = 0 + 0,25 = 0,25
у1 = у0 + 0,5h [(2x0 – y20 + х40) + (2x1 – ( y1*)2 + х41)] =
= 0 + 0,5·0,25·[(2∙0 – 02 + 04) + (2∙0,25 – 02 + 0,254) = 0,063;
Второй шаг: i = 1, x1 = 0,25, y1 = 0,063
23
y*2= y1 + h (2x1 – y21 + х41) = 0,063 + 0,25·(2∙0,25 – 0,0632 + 0,254) = 0,188;
x2 = x1 + h = 0,25 + 0,25 = 0,5
у2 = у1 + 0,5h [(2x1 – y21 + х41) + (2x2 – ( y2*)2 + х42)] =
= 0,063 + 0,5·0,25·[(2∙0,25 – 0,0632 + 0,254) + (2∙0,5 – 0,1882 + 0,54)] =
0,254;
Третий шаг: i = 2, x2 = 0,5, y2 = 0,254
y*3= y2 + h (2x2 – y22 + х42) = 0,254 + 0,25·(2∙0,5 – 0,2542 + 0,54) = 0,5035:
x3 = x2 + h = 0,5 + 0,25 = 0,75
у3 = у2 + 0,5h f [(2x2 – y22 + х4 2) + (2x3 – ( y3*)2 + х43)] =
= 0,254 + 0,5·0,25·[(2∙0,5 – 0,2542 + 0,54) + (2∙0,75 – 0,50352 + 0,754)] =
0,574.
7.3.3 Усовершенствованный метод Эйлера
Сущность усовершенствованного метода Эйлера состоит в следующем:
сначала вычисляются вспомогательные значения искомой функции yi+1/2 в
точках xi+1/2 = xi + h/2 помощью формулы
y
i1/2
h
 y  y' ,
i
2 i
'
затем находят значение f (x,y) в средней точке yi1/2
 f (xi1/2 , yi1/2 ) и
определяют
y
i1
 y  h  y'
i
.
i1/2
Оценка погрешности в точке хi может быть получена с помощью
«двойного просчета»: расчет повторяют с шагом h/2 и погрешность более
точного y*i ( при шаге h/2) оценивают приближенно по формуле
24
y*  y(x ) 
i
i
1
3
yi*  y i ,
где у(х) – точное решение ДУ.
Пример 7.5.
Проинтегрировать усовершенствованным методом Эйлера ДУ y’ = y – x
при начальных условиях х0 = 0, у0 = 1,5 на отрезке [0, 1], приняв h = 0,25.
Первый шаг: i = 0, x0 = 0, y0 = 1,5.
y’0+1 = y0 - x0 = 1,5 – 0 = 1,5;
x0+1/2 = x0 +
h
0
2
0, 25
h
y '
21
 0,125; y0+1/2 = y0 +
2
h
2
0, 25
1,5  0,1875;
2
y '  1,5  0,1875  1,6875 ;
1
'
y’0+1/2 = y0+1/2 – x0+1/2 = 1,6875 – 0,125 = 1,5625; h  y
01/ 2
 0,25 1,5625 
0,3906;

 y  h  y'
y
01
0
 1,5 + 0,3906 = 1,8906;
01/2
Второй шаг: i = 1, x1 = x0 + h = 0 + 0,25 = 0,25; y1 = 1,8906.
y’2 = y1 – x1 = 1,8906 – 0,25 = 1,6406;
h
2
h
x1+1/2 = x1 +
y1+1/2 = y1 +
h
2
 0, 25 
0, 25
y '
2
0, 25
1,6406  0,2051;
2
 0,3755 ;
2
2
'
y  1,8906  0, 2051  2,0957 ;
2
'
y’1+1/2 = y1+1/2 – x1+1/2 = 2,0957 – 0,375 = 1,7207; h  y
11/ 2
 0,25 1,7207 
0,4302;

y  y  h  y'
2
1
 1,8906 + 0,4302 = 2,3208.
11/2
25
За повышение точности приходится расплачиваться дополнительными
*
затратами времени на вычисление функции yi1
.
Более высокая точность может быть получена, если улучшить
аппроксимацию производной, сохраняя большее число членов ряда Тейлора.
7.3.4 Методы Рунге-Кутта
Метод Рунге-Кутта является одним из методов повышенной точности.
Он имеет много общего с методом Эйлера. Метод Эйлера можно считать
методом Рунге-Кутта первого порядка ( в разложении в ряд Тейлора остается
только первая производная).
Пусть требуется найти численное решение уравнения y’ = f (x, y) на
отрезке [a, b] с начальными условиями у(х0) = у0.
Разобьем отрезок [a, b] на n равных частей с точками xi = x0 + i·h (i = 0,
1,…, n), где h = (b - a)/n – шаг интегрирования. В методе Рунге-Кутта, так же
как и в методе Эйлера, последовательные значения yi искомой функции у
определяются по формуле
yi +1 = yi + Δ yi.
Если разложить функцию у в ряд Тейлора и ограничиться членами до h4
включительно, то приращение функции Δ y можно представить в виде
y  y(x  h)  y(x)  hy' (x) 
h2
2
y (x) 
''
h3
6
y (x) 
'''
h4
24
y'''' (x) ,
(7.2)
где производные y'' (x), y''' (x), y'''' (x) находят последовательным
дифференцированием из уравнения y’ = f (x, y).
Вместо непосредственных вычислений по формуле (7.2) в методе РунгеКутта определяют четыре числа:
k1  hf (x, y);
h
k
k  hf (x  , y  1 );
2
2
2
h
k2
k  hf (x  , y  );
3
2
2
k4  hf (x  h, y  k3 ).
(7.3)
26
Можно доказать, что если числам k1, k2, k3, k4 придать соответственно
вес 1/6; 1/3; 1/3; 1/6, то средневзвешенное этих чисел, т.е.
1
1
1
1
k k  k  k
6 1 3 2 3 3 6 4
с точностью до четвертых степеней равно значению Δу, вычисленному по
формуле (7.2)
1
y  (k 2k 2k k ).
2
3
4
6 1
(7.4)
Таким образом, для каждой пары текущих значений xi и yi по формулам
(7.3) определяют значения
ki  hf (x , y );
1
i
i

i
h
1k
k  hf (xi  , yi  );
2
2
i
2

(7.5)

h
ki
k3i  hf (xi  , y  2 );
2
2
ki  hf (x  h, y  ki ).
4
i
i
3
по формуле (7.4) находят
1
y  (ki  2ki  2ki  ki ).
i
2
3
4
6 1
Метод Рунге-Кутта имеет порядок точности h4 на всем отрезке [a, b].
Оценка точности этого метода затруднительна. Грубую оценку погрешности
можно получить с помощью «двойного просчета» по формуле
y i  y(xi ) 
*
y*  y
i
i
15
,
27
где у(хi) – значение точного решения ДУ в точке хi, а yi* и уi – приближенные
значения, полученные с шагом h/2 и h.
Если ε – заданная точность решения, то число n (число делений) для
определения шага интегрирования h = (b - a)/n выбирается таким образом,
чтобы
h4 < ε.
Однако шаг расчета можно менять при переходе от одной точки к
другой. Для оценки правильности выбора шага h используется равенство
k2i k 3i
q i
,
k1  k2i
где q должно быть равно нескольким сотым, в противном случае шаг h
уменьшают.
Пример 7.6.
Дано ДУ у’ = у – х, удовлетворяющее начальному условию у(0) = 1,5.
Вычислить с точностью до 0,01 решение этого уравнения на интервале [0,
1,5]. Вычисление выполнить методом Рунге-Кутта с двумя запасными
знаками.
Выбираем начальный шаг из условия h4 < ε = 0,01:
h  4 0,01  0,316.
Для удобства вычислений примем h = 0,25.
Первый шаг: i = 0, x0 = 0, y0 = 1,5.
k0  h( y  x )  0,25(1,5  0)  0,375;
1
0
0
0,25 )]  0,3906;
1
k)0  (x  h)]  0, 25[(1,5  0,375

h[(
y

)  (0 
0
0
k
2
2
2
2
0
2
0
h
2k

h[(
y

)
 (x0  )]  0, 25[(1,5  0,3906)  (0  0,25 )]  0,3926;
0
k
2
2
2
2
0
3
28
k 0  h[( y  k 0 )  (x  h)]  0,25[(1,5  0,3926)  (0  0,25)]  0,4419;
4
y 
0
0
3
0
1
1
(k 0  2k 0  2k 0  k 0 )  (0,375  2  0,3906  2  0,3926  0, 4419)  0,3972;
2
3
4
6 1
6
y1  y0  y0 1,5  0,3972 1,8972 .
Второй шаг: i = 1, x1 = x0 + h = 0 + 0,25 = 0,25, y1 = 1,8972.
k1  h( y  x )  0,25(1,8972  0,25)  0,4118;
1
1
1
k11
h
0, 4118
0, 25
k2  h[( y1  )  (x1  )]  0, 25[(1,8972 
)  (0, 25 
)]  0, 4320;
2
2
2
2
1
k21
h
0, 4320
0, 25
k3  h[( y1  )  (x1  )]  0, 25[(1,8972 
)  (0, 25 
)]  0, 4346;
2
2
2
2
1
k1  h[( y  k1)  (x  h)]  0,25[(1,8972  0,4346)  (0,25  0,25)]  0,4580;
4
0
3
1
1
1
y  (k1  2k1  2k1  k1 )  (0, 4118  2  0, 4320  2  0, 4346  0, 4580)  0, 4338;
2
2
3
4
6 1
6
y2  y1  y2 1,8972  0,4338  2,331.
Метод Рунге-Кутта может быть применен к решению систем ДУ.
Пусть дана система ДУ первого порядка:
 y '  f1 (x, y, z),
 '
 z  f2 (x, y, z)
с начальными условиями х = х0, у(х0) = у0, z(х0) = z0.
В этом случае параллельно определяются числа Δуi, Δzi:
1
y  (ki  2ki  2ki  ki ),
i
2
3
4
6 1
29
1
z  (li  2li  2li  li ),
i
2
3
4
6 1
где
k i  hf (x , y , z );
1
1
i
i
i
l1i  hf2(x , y , z );
i
i
i
ki 1 li
k2  hf1(xi  , yi  , zi  );
2
2
2
h
i
1
ki 1 li
l2  hf2 (xi  , yi  , zi  );
2
2
2
h
i
1
i
i
k
l
k3  hf1(xi  , yi  2 , zi  2 );
2
2
2
h
i
ki 2 li
l3  hf2 (xi  , yi  , zi  );
2
2
2
h
i
2
ki  hf (x  h, y  ki , z  li );
4
1
i
i
3
i
3
li  hf (x  h, y  ki , z  li ).
4
2
i
i
3
i
3
Тогда получим решение системы
yi1  yi  yi ,
zi1  zi  zi .
Пример 7.7.
Задана система ДУ
 ' 2yx
y

z

z'  2 y
zx

при начальных условиях х0 = 0,5; у0 = 1,0; z0 = 1,0. Найти решение
системы при х = 0,6. Вычисления вести с пятью знаками после запятой.
30
Выбираем шаг h = 0,1 и определим коэффициенты:
k1 h
2 y0  x0  0,12 1 0,5  0,15000;
1,0
z0
l  h 2 y0  0,1 2 1  0,13333;
z x
1,0  0,5
1
0
2( y 
k1
h
)  (x  )
2
0
0
k2  h[
0
2(1 
2
)  (0,5 
0.1
2)
2
]  0,1[
z0  l1
2
0,15000
]  0,15000;
1,0 0,13333
2
0,15000
2( y0  k1 )
2(1

)
2
l2 h[
l1
h ]  0,1[
0,13333 2
0,1 ]  0,13299;
(z  )  (x  )
(1,0 
)  (0,5 
0
2( y 
k2
0
k3  h[
)
0
2
2
h
)  (x  )
2
0
2
2(1 
)  (0,5 
2
]  0,1[
z0  2l
2
2
0,15000
2
0.1
2 )
]  0,15002;
1,0 0,13299
2
0,15000
2( y0  k2 )
2(1  2 )
2
l  h[
l
h ]  0,1[
0,13299
0,1 ]  0,13300;
(z  2 )  (x  )
(1,0 
)  (0,5 
3
0
k4  h[
0
2
2
2
2( y0  k3 )  (x0  h)
2(1 0,15002)  (0,5  0,1)
]

0,1[
]  0,15005;
z l
1,0  0,13300
0
l3 h[
)
2
3
2(1 0,15002)
2( y0  k3 )
]  0,1[
]  0,13272;
(z  l )  (x  h)
(1,0  0,13300)  (0,5  0,1)
0
3
0
Вычисляем приращение искомых функций
1
1
y  (k  2k  2k  k )  (0,15  2  0,15000  2  0,15002  0,15005) 
0
2
3
4
6 1
6


31
 0,150015,
32
1
1
z  (l  2l  2l  l )  (0,13333  2  0,13299  2  0,13300  0,13272) 
0
2
3
4
6 1
6
 0,133005.


Значение искомых функций в точке х = 0,6
y1  y0  y0 1 0,150015 1,150015,
z1  z0  z0 1 0,133005 1,133005.
7.3.5 Общая характеристика одношаговых методов
Одношаговым методам присущи следующие общие черты:
1.
Чтобы получить информацию в новой точке надо иметь данные
лишь в одной предыдущей точке. Это свойство называется
«самостартованием» - «self-starting behavior».
2.
В основе всех одношаговых методов лежит разложение функции
в ряд Тейлора, в котором сохраняются члены, содержащие шаг h в степени до
m включительно. Целое число m называется порядком метода.
3.
Все одношаговые методы не требуют действительного
вычисления производных – вычисляется лишь сама функция, однако могут
потребоваться значения в нескольких промежуточных точках.
4.
Свойство «самостартования» позволяет легко менять величину
шага.
Документ
Категория
Без категории
Просмотров
3
Размер файла
188 Кб
Теги
ssu, numerical, system
1/--страниц
Пожаловаться на содержимое документа