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

Содержание

Численные методы решения обыкновенных дифференциальных уравнений

Тема 7

Численные методы решения

обыкновенных дифференциальных уравнений

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

Основы теории обыкновенных

дифференциальных уравнений

Обыкновенным дифференциальным уравнением называется уравнение вида

                       (3.1)

связывающее одну независимую переменную , искомую функцию  и ее производные вплоть до n-го порядка. Порядком обыкновенного дифференциального уравнения называется порядок старшей производной от искомой функции.

Дифференциальное уравнение называется линейным, если оно имеет вид

Так, например,  – линейное дифференциальное  уравнение  второго порядка;  – нелинейное уравнение.

Общим интегралом уравнения (3.1) называют функцию

               (3.2)

связывающую независимую переменную , искомую функцию  и п постоянных интегрирования  с помощью уравнения . Таким образом, функция  входит в функцию (3.2) неявным образом, причем число постоянных интегрирования равно порядку уравнения.

Общим решением обыкновенного дифференциального уравнения называется функция

              (3.3)

связывающая независимую переменную и п постоянных интегрирования, иначе говоря, уравнение (3.3) определяет функцию   явным образом.

Для определения постоянных интегрирования задаются дополнительные условия. Их число равно числу постоянных интегрирования. Если в дополнительные условия подставить функцию (3.2) и решить полученную систему относительно постоянных интегрирования  а затем полученные решения подставить в (3.2), то получим частный интеграл .

Аналогичные процедуры с общим решением (3.3) дают частное решение

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

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

Таким образом, если задать числа  то задача Коши для уравнения (3.1) имеет вид

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

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

Во многих практических приложениях независимая переменная обозначается через  и имеет смысл времени, задача Коши называется «начальной».

Актуальность задачи Коши для многих областей науки и техники явилась причиной разработки для ее решения большого количества методов. Рассматриваются две группы численных методов решения задачи Коши.

– Одношаговые методы, в которых для нахождения решения в некоторой точке отрезка используется информация лишь в одной предыдущей точке. К этим методам относятся методы Эйлера и Рунге – Кутты.

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

К таким методам относится метод Адамса.

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

  тогда,

Пример 

Записать в нормальной форме Коши

Обозначим:

Понятие о методе конечных разностей.

Порядок разностной схемы

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

– Область непрерывного изменения аргумента (или аргументов) заменяется конечным дискретным множеством точек, называемых разностной сеткой. В разностной сетке выделяются внутренние и граничные узлы. Решение разыскивается во внутренних узлах, а в граничных узлах значение искомой функции задается при аппроксимации граничных условий исходной дифференциальной задачи. Функция дискретного аргумента, определенная на разностной сетке, называется сеточной функцией.

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

– Осуществляется решение системы алгебраических уравнений каким-либо из известных методов.

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

Методы сеток делятся на явные и неявные:

– Явные – когда значение  на  шаге определяется явно:

, где -некоторая функция, зависящая от конкретного  метода (кроме последней рассчитанной точки ) могут использоваться еще  предыдущих значений;

– Неявные – когда искомая величина  входит как в левую, так и в правую часть:

 

Явные и неявные методы делятся на 1-шаговые и многошаговые:

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

В шаговых – требуется информация о предыдущих точках.

В случае нелинейных систем итерационные процедуры, как правило, сводят их к линейным системам.

Продемонстрируем эти этапы на примере решения задачи Коши для простейшего линейного уравнения первого порядка. Пусть требуется найти решение задачи Коши на интервале (0, 1) для уравнения

             (3.4)

Непрерывный интервал можно заменить множеством точек . Расстояние между точками  называется шагом разностной сетки. В общем случае эти шаги могут быть переменными. Рассмотрим случай постоянного шага и примем . Заменим дифференциальное уравнение разностным в точке :

                                                             (3.5)

Таким образом, первая производная аппроксимирована односторонней разностью вперед с первым порядком точности. Перепишем уравнение (3.5) в виде

                                                                                               (3.6)

Здесь  Из (3.6) следует, что

.                (3.7)

Из рекуррентного соотношения (3.7) имеем точное решение разностного уравнения (3.6)

             (3. 8)

Очевидно, что точное решение разностных уравнений можно получить лишь для простейших случаев. Точное решение задачи (3.4) имеет вид

.                (3.9)

Определим погрешность разностного решения. Имеем

. (3.10)

Предполагая шаг разностной схемы малым, представим член  в виде

 

 Подставляя это выражение в (3.10), получим

.                                                                          (3.11)

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

Говорят, что разностная схема (3.5) имеет первый порядок точности на интервале. Можно показать, что на шаге эта разностная схема имеет второй порядок точности. Представим точное решение в виде ряда Тейлора в точке. Имеем

.                                                                           (3.12)

Поскольку

,

то

.                                                                            (3.13)

Из соотношений (3.7) и (3.13) следует, что погрешность на шаге имеет порядок . Очевидно, что погрешность на интервале выше на порядок, так как она накапливается по мере увеличения числа шагов. Отметим еще, что погрешность решения стремится к нулю при , т. е. решение разностной задачи сходится к точному. Однако для сходимости необходимы не только аппроксимация дифференциального уравнения, но и устойчивость разностной схемы.

Численный метод называется устойчивым, если численные результаты непрерывно зависят от входных данных и если погрешность остается ограниченной при заданных приделах изменения параметров численного алгоритма (шагов сетки, числа итераций).

Сходимость является основной проблемой, от которой зависит точность всей задачи.

Численный алгоритм называется сходящимся, если при стремлении его параметров к предельным значениям, например при , (итерационный шаг при – число итераций), результаты стремятся к точному решению.

Метод Эйлера.

Метод Эйлера с пересчетом

Метод Эйлера является простейшим методом решения задачи Коши и имеет невысокую точность, поэтому на практике его используют достаточно редко. Однако на его основе в дальнейшем легче объяснить алгоритмы более эффективных (и, как правило, более сложных) методов и способы построения и исследования этих алгоритмов.

Пусть дано дифференциальное уравнение

                   (1)

с начальным условием

.

Выбрав достаточно малый шаг , построим систему равноотстоящих точек

.                                                                                      (2)

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

с вершинами , звенья которой  прямолинейны между прямыми и имеют подъем

 

                        (3)

(так называемая ломаная Эйлера).

Таким образом, звенья  ломаной Эйлера в каждой вершине имеют направление, совпадающее с направлением уд интегральной кривой уравнения (1), проходящей через точку.

Из формулы (3) вытекает, что значения  могут быть определены (метод Эйлера) по формулам

.                                                                         (4)

Для геометрического построения ломаной Эйлера выберем полюс  и на оси ординат отложим отрезок  (рис. 28). Очевидно, угловой коэффициент луча  будет равен ; поэтому, чтобы получить первое звено ломаной Эйлера, достаточно из точки  , провести прямую , параллельную лучу, до пересечения с прямой  в некоторой точке .

Приняв точку  за исходную, откладываем на оси ординат отрезок  и через точку проводим прямую до пересечения в точке с прямой  и т. д.

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

  1.  малая точность;
  2.  систематическое накопление ошибок.

Можно доказать, что если правая часть  уравнения (1) непрерывна, то последовательность ломаных Эйлера при  на достаточно малом отрезкеравномерно стремится к искомой интегральной кривой.

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

Пример. Применяя метод Эйлера, составить на отрезке [0,1] таблицу значений интеграла дифференциального уравнения

                                                  ,                                             (5)                                                          

удовлетворяющего начальному условию , выбрав шаг .

Решение. Результаты вычислений приведены в таблице 35. Для сравнения в последнем столбце помещены значения точного решения

.

Из приведенной таблицы видно, что абсолютная погрешность значения  составляет . Отсюда относительная погрешность примерно равна 3%.

Для сравнения приводим график точного решения (выделенный жирной линией) и соответствующую ломаную Эйлера  (рис. 29).

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

,

т.е. для этого отрезка имеется погрешность порядка .

Таблица 35

Интегрирование дифференциального уравнения (5) методом Эйлера

Точное значение

0

1

2

3

4

5

6

7

8

9

10

 0

0,1

0,2

0,3

0,4

0,5

0,6

0,7

0,8

0,9

1,0

  1

  1

  1,005

1,0151

1,0303

1,0509

1,0772

1,1095

1,1483

1,1942

1,2479

        0

        0,5

0,1005

0,1523

0,2067

0,2627

0,3232

0,3883

0,4593

0,5374

        0

        0,005

0,0101

0,0152

0,0206

0,0263

0,0323

0,0388

0,0459

0,0537

        1

1,0025

1,0100

1,0227

1,0408

1,0645

1,0942

1,1303

1,1735

1,2244

1,2840

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

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

Подобные схемы часто называют схемами типа прогноз – корректор (предиктор – корректор). Сначала определяется «грубое приближение» решения

,

исходя из которого находится направление поля интегральных кривых

.

 Затем решение уточняют, приближенно полагают

                                                                                              (5)

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

,

построим итерационный процесс

                                 .                      (7)

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

,

где –  общая часть приближений  и .

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

Отметим в заключение, что метод Эйлера  с итерационной обработкой ординат дает на каждом шаге погрешность порядка и нередко применяется в вычислительной практике.

Пример 2. Применяя метод итерационной обработки, с точностью до четырех совпадающих десятичных знаков найти значение

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

.

Решение. Примем шаг

.

Применяя итерационный процесс (7) и учитывая, что

,

последовательно будем иметь

Следовательно, удержав запасной знак, можно положить

.

Аналогично, приняв  и  за исходные данные и принимая во внимание, что

,

С помощью того же итерационного процесса (7) получим

Отсюда

.

Для сравнения приводим точное значение

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

Эйлера с последующей итерацией значений

Точное значение

0

1

2

3

4

5

6

7

8

9

10

 0

0,1

0,2

0,3

0,4

0,5

0,6

0,7

0,8

0,9

1,0

  1

  1,0025

  1,0100

1,0227

1,0408

1,0646

1,0943

1,1305

1,1738

1,2248

1,2812

        0

        0,0501

0,1010

0,1534

0,2082

0,2661

0,3283

0,3957

0,4695

0,5512

        0

        0,0050

0,0101

0,0153

0,0208

0,0266

0,0328

0,0396

0,0470

0,0551

        1

1,0025

1,0100

1,0227

1,0408

1,0645

1,0942

1,1303

1,1735

1,2244

1,2840

Метод Рунге – Кутты

 Метод Рунге — Кутты позволяет строить схемы различного порядка точности. Основная идея метода состоит в построении специального алгоритма – такого, чтобы приращение функции на шаге  совпадало с приращением , которое определяется из ряда Тейлора (3.15) с учетом возможно большего числа членов. При этом вторые и следующие производные определяются не дифференцированием, а путем многократного вычисления функции  в некоторых промежуточных точках между и .

Пусть дано дифференциальное уравнение первого порядка

                                                                                                                         (1)

с начальным условием

 

 

 Выберем шаг  и для кратности введем обозначения  и .

 

Рассмотрим числа:

                                                                                (2)

Согласно обычному методу Рунге – Кутта, последовательные значения   искомой функции  определяются по формуле

,

где

                                        .                      (3)

Формула (9) имеет четвертый порядок точности. Получены также формулы типа Рунге – Кутта с иными порядками точности.

Для вычисления по формуле (3) удобно пользоваться схемой, приведенной в таблице 39.

Эффективная оценка погрешности метода Рунге – Кутта затруднительна. Поэтому для определения правильности выбора шага h на практике обычно на каждом этапе из двух шагов применяют двойной пересчет. А именно, исходя из текущего верного значения ,  вычисляют величину  двумя способами: один раз с шагом h, а другой раз с двойным шагом . Если расхождение полученных значений не превышает допустимой погрешности, то шаг h для данного этапа выбран правильно и полученное  с  его помощью значение можно принять за . В противном случае шаг уменьшаю в два раза.

Схема метода Рунге – Кутта

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

Пример 1. Методом Рунге – Кутта вычислить на отрезке [0; 0,5] интеграл дифференциального уравнения

                                                                                               (10)

приняв шаг h = 0,1.

Решение. Покажем начало процесса.

Вычисление .

Последовательно имеем

Отсюда

 

и, следовательно,

                      .

Аналогично вычисляются дальнейшие приближения. Результаты вычислений приведены в таблице 40.

Таким образом, у (0,5) = 1,7974.

Для сравнения приводим точное решение:

,

откуда

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

Таблица 40

Интегрирование дифференциального уравнения (10)

методом Рунге – Кутта

0

           0

0,05

0,05

           0,1

        1

        1,05

        1,055

1,1105

         0,1

         0,11

0,1105

0,1210

1

           0,1

0,15

0,15

           0,2

        1,1103

        1,1708

        1,1763

1,2429

         0,1210

         0,1321

0,1326

0,1443

2

           0,2

0,25

0,25

           0,3

        1,2427

        1,3149

        1,3209

1,3998

         0,1443

         0,1565

0,1571

0,1700

3

           0,3

0,35

0,35

           0,4

        1,3996

        1,4846

        1,4904

1,5836

         0,1700

         0,1835

0,1840

0,1984

4

           0,4

0,45

0,45

           0,5

        1,5836

        1,6828

        1,6902

1,7976

         0,1984

         0,2133

0,2140

0,2298

5

0,5

1,7974

Многошаговый метод Адамса

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

Общая схема построения многошаговых методов выглядит следующим образом. Пусть нам известно приближенное решение.

Изложим метод Адамса применительно к уравнению первого порядка

                                                                                                             (1)

с начальным условием

                                                                                                                 (2)

Пусть  – система равноотстоящих значений с шагом h и . Очевидно, имеем

 . (3)

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

, (4)

где

,

или

                          , (4’)

Подставляя выражение (4′) в формулу (3) и учитывая, что

будем иметь

.

Отсюда получаем экстраполяционную формулу Адамса

. (5)

Для начала процесса нужны четыре начальных значения:

так называемый начальный отрезок, который определяется, исходя из начального условия (2), каким–нибудь  численным методом.

Можно, например, использовать метод Рунге – Кутта  или разложение в ряд Тейлора

.,

где (или с соответствующим изменением нумерации). Зная эти значения, из уравнения (1) можно найти значения производных

и составить таблицу разностей:

                                                      (6)

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

Для работы на электронных счетных машинах формулу Адамса (5) выгодно применять в раскрытом виде. Учитывая, что

после приведения подобных членов имеем

,

причем

Таким образом, метод Адамса имеет четвертый порядок точности на интервале. Чтобы начать счет по схеме Адамса, необходимо знать решение в четырех начальных точках .

По существу, интерполяционный многочлен  в формуле (3.32) используется вне области интерполяции, т. е. в данном случае это экстраполяционный многочлен. Однако, поскольку интервал  мал, ошибка за счет экстраполяции невелика. Недостающие значения функции вычисляются в точках , как правило, по методу Рунге – Кутты  соответствующего порядка. Это является недостатком метода, так как увеличивает объем программы для компьютера. Преимущество многошаговых методов заключается в том, что на каждом шаге правая часть дифференциального уравнения вычисляется только один раз, а в методе Рунге – Кутты четвертого порядка точности на каждом шаге функция  вычисляется четыре раза.

Достоинство метода Адамса по сравнению с методом Рунге – Кутты  заключается в простоте оценки остаточного члена метода.

Метод Адамса распространяется и на системы дифференциальных уравнений 1-го порядка.

Пример 1. Методом Адамса найти на отрезке [0,1] интеграл уравнении

.

Решение. Примем шаг . Для начала процесса используем значения, найденные методом Рунге – Кутта (см. §7, пример 1), т. е.

.

Дальнейшие вычисления располагаем в двух бланках: основном (таблица 42) и вспомогательном (таблица 43).

Если правая часть дифференциального уравнения сложна, то в основном бланке приходится вводить промежуточные графы. При заполнении вспомогательного бланка используются согласно формуле (5) наклонные строки.

Таблица  42

Основной бланк для интегрирования дифференциального уравнения

методом Адамса

0

1

2

3

4

5

6

7

8

9

10

0

0,1

0,2

0,3

0,4

0,5

0,6

0,7

0,8

0,9

1,0

 1

1,1103

1,2427

1,3996

1,5834

1,7971

2,0440

2,3273

2,6508

3,0190

3,4362

0,1838

0,2137

0,2469

0,2833

0,3235

0,3682

0,4172

0,1000

0,1210

0,1443

0,1700

0,1983

0,2297

0,2644

0,3027

0,3451

0,3919

210

233

257

283

314

347

383

424

468

23

24

26

31

33

36

41

44

1

2

5

2

3

5

3

 1

1,1103

1,2428

1,3997

1,5836

1,7974

2,0442

2,3275

2,6511

3,0192

3,4366

Вспомогательный бланк для интегрирования дифференциального

уравнения методом Адамса

3

4

5

6

7

8

9

0,1700

    128

     10

       0

0,1983

   142

     11

      1

0,2297

   157

    13

      2

0,2644

    174

     14

       1

0,3027

   192

     15

      1

0,3451

   212

    17

     2

0,3919

    234

      18

       1

0,1838

0,2137

0,2469

0,2833

0,3235

0,3682

0,4172

В последнем столбце приведены для сравнения точные значения решения

.

Отсюда видно, что максимальная ошибка приближенного решения у не превосходит четырех единиц последнего десятичного разряда.

Неявные схемы

Рассмотрим следующую разностную схему для численного решения уравнения (3.14) при условии, что вычислено значение функции u в точке . Имеем

           ,                       (3.34)

где s — параметр разностной схемы такой, что . Говорят, что разностная схема является явной, если правая часть дифференциального уравнения записывается в известной точке . Схема называется неявной, если правая часть дифференциального уравнения записывается в неизвестной точке . Схема называется полуявной (полунеявной), если правая часть с определенными весами записывается в известной и неизвестной точках. В соответствии с этими терминами схема (3.34) является явной при s = 1, неявной при s = 0 и полуявной (полунеявной) при . В явных схемах решение в точке  определяется непосредственно из алгебраических соотношений с известными коэффициентами. Все рассмотренные выше схемы методов Эйлера, Рунге – Кутты, Адамса (кроме схемы метода Эйлера с пересчетом) являются явными. В неявных и полунеявных схемах при  для нахождения решения в точке  необходимо решить в общем случае нелинейное уравнение (3.34).

Нелинейное уравнение (3.34) можно разрешить различными методами, например, методами простой итерации, дихотомии, Ньютона и т. д. В рассматриваемом случае логично использовать метод простой итерации, поскольку нелинейное уравнение (3.34) уже разрешено относительно. Тогда, если номер итерации обозначить верхним индексом, итерационный процесс можно записать в виде

              .                        (3.35)

Известно, что сходимость метода простой итерации зависит от выбора начального приближения. В рассматриваемом случае в качестве   можно выбрать значение  или какое-либо другое, например, значение, вычисленное по методу Эйлера, достаточно близкое к искомому значению . Вычисления по формуле (3.35) продолжаются до получения заданной точности .

Таким образом, очевидно, что при расчете на шаге явные схемы требуют значительно меньшего числа операций, чем неявные. Однако неявные схемы значительно более устойчивы, чем явные. Поэтому в неявных схемах ограничения на шаг разностной схемы,  связанные с устойчивостью, как правило, менее жесткие, чем в явных. В результате может оказаться, что общее число операции при заданной точности и при расчетах на интервале интегрирования в неявных схемах меньше чем в явных.

Pers.narod.ru. Обучение. Лекции по численным методам. Численное решение обыкновенных дифференциальных уравнений

Pers.narod.ru. Обучение. Лекции по численным методам. Численное решение обыкновенных дифференциальных уравнений

Этот сайт больше не обновляется. Подключите Javascript, чтобы увидеть новый адрес страницы или перейдите к статье

5. Численное решение обыкновенных дифференциальных уравнений

Численные методы решения задачи Коши для ОДУ первого порядка
Численные методы решения систем ОДУ первого порядка
Метод конечных разностей решения краевых задач для ОДУ

Многие задачи науки и техники сводятся к решению обыкновенных дифференциальных уравнений (ОДУ). ОДУ называются такие уравнения, которые содержат одну или несколько производных от искомой функции. В общем виде ОДУ можно записать следующим образом:

, где x – независимая переменная, – i-ая производная от искомой функции. n - порядок уравнения. Общее решение ОДУ n–го порядка содержит n произвольных постоянных , т.е. общее решение имеет вид .

Для выделения единственного решения необходимо задать n дополнительных условий. В зависимости от способа задания дополнительных условий существуют два различных типа задач: задача Коши и краевая задача. Если дополнительные условия задаются в одной точке, то такая задача называется задачей Коши. Дополнительные условия в задаче Коши называются начальными условиями. Если же дополнительные условия задаются в более чем одной точке, т.е. при различных значениях независимой переменной, то такая задача называется краевой. Сами дополнительные условия называются краевыми или граничными.

Ясно, что при n=1 можно говорить только о задачи Коши.

Примеры постановки задачи Коши:

Примеры краевых задач:

Решить такие задачи аналитически удается лишь для некоторых специальных типов уравнений.

Численные методы решения задачи Коши для ОДУ первого порядка

Постановка задачи. Найти решение ОДУ первого порядка

на отрезке при условии

При нахождении приближенного решения будем считать, что вычисления проводятся с расчетным шагом , расчетными узлами служат точки промежутка [x0, xn].

Целью является построение таблицы

xi

x0

x1

xn

yi

y0

y1

yn

т.е. ищутся приближенные значения y в узлах сетки.

Интегрируя уравнение на отрезке , получим

Вполне естественным (но не единственным) путем получения численного решения является замена в нем интеграла какой–либо квадратурной формулой численного интегрирования. Если воспользоваться простейшей формулой левых прямоугольников первого порядка

,

то получим явную формулу Эйлера:

, .

Порядок расчетов:

Зная , находим , затем т.д.

Геометрическая интерпретация метода Эйлера:

Пользуясь тем, что в точке x0 известно решение y(x0) = y0 и значение его производной , можно записать уравнение касательной к графику искомой функции в точке : . При достаточно малом шаге h ордината этой касательной, полученная подстановкой в правую часть значения , должна мало отличаться от ординаты y(x1) решения y(x) задачи Коши. Следовательно, точка пересечения касательной с прямой x = x1 может быть приближенно принята за новую начальную точку. Через эту точку снова проведем прямую , которая приближенно отражает поведение касательной к в точке . Подставляя сюда (т.е. пересечение с прямой x = x2), получим приближенное значение y(x) в точке x2: и т.д. В итоге для i–й точки получим формулу Эйлера.

Явный метод Эйлера имеет первый порядок точности или аппроксимации.

Если использовать формулу правых прямоугольников: , то придем к методу

, .

Этот метод называют неявным методом Эйлера, поскольку для вычисления неизвестного значения по известному значению требуется решать уравнение, в общем случае нелинейное.

Неявный метод Эйлера имеет первый порядок точности или аппроксимации.

Модифицированный метод Эйлера: в данном методе вычисление состоит из двух этапов:

Данная схема называется еще методом предиктор – корректор (предсказывающее – исправляющее). На первом этапе приближенное значение предсказывается с невысокой точностью (h), а на втором этапе это предсказание исправляется, так что результирующее значение имеет второй порядок точности.

Методы Рунге – Кутта: идея построения явных методов Рунге–Кутты p–го порядка заключается в получении приближений к значениям y(xi+1) по формуле вида

,

где

…………………………………………….

.

Здесь an, bnj, pn, – некоторые фиксированные числа (параметры).

При построения методов Рунге–Кутты параметры функции (an, bnj, pn) подбирают таким образом, чтобы получить нужный порядок аппроксимации.

Схема Рунге – Кутта четвертого порядка точности:

Пример. Решить задачу Коши:

.

Рассмотреть три метода: явный метод Эйлера, модифицированный метод Эйлера, метод Рунге – Кутта.

Точное решение:

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

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

Расчетные формулы метода Рунге – Кутта:

x

y1

y2

y3

точное

0

1.0000

1.0000

1.0000

1.0000

0.1

1.2000

1.2210

1.2221

1.2221

0.2

1.4420

1.4923

1.4977

1.4977

0.3

1.7384

1.8284

1.8432

1.8432

0.4

2.1041

2.2466

2.2783

2.2783

0.5

2.5569

2.7680

2.8274

2.8274

0.6

3.1183

3.4176

3.5201

3.5202

0.7

3.8139

4.2257

4.3927

4.3928

0.8

4.6747

5.2288

5.4894

5.4895

0.9

5.7377

6.4704

6.8643

6.8645

1

7.0472

8.0032

8.5834

8.5836

y1 – метод Эйлера, y2 – модифицированный метод Эйлера, y3 – метод Рунге Кутта.

Видно, что самым точным является метод Рунге – Кутта.

Численные методы решения систем ОДУ первого порядка

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

Покажем это для случая системы двух уравнений первого порядка:

Явный метод Эйлера:

Модифицированный метод Эйлера:

Схема Рунге – Кутта четвертого порядка точности:

К решению систем уравнений ОДУ сводятся также задачи Коши для уравнений высших порядков. Например, рассмотрим задачу Коши для уравнения второго порядка

Введем вторую неизвестную функцию . Тогда задача Коши заменяется следующей:

Т.е. в терминах предыдущей задачи: .

Пример. Найти решение задачи Коши:

на отрезке [0,1].

Точное решение:

Действительно:

Решим задачу явным методом Эйлера, модифицированным методом Эйлера и Рунге – Кутта с шагом h=0.2.

Введем функцию .

Тогда получим следующую задачу Коши для системы двух ОДУ первого порядка:

Явный метод Эйлера:

Модифицированный метод Эйлера:

Метод Рунге – Кутта:

Схема Эйлера:

X

y

z

y теор

z теор

y-y теор

0

1

0

1

0

0

0.2

1

-0.2

0.983685

-0.14622

0.016315

0.4

0.96

-0.28

0.947216

-0.20658

0.012784

0.6

0.904

-0.28

0.905009

-0.20739

0.001009

0.8

0.848

-0.2288

0.866913

-0.16826

0.018913

1

0.80224

-0.14688

0.839397

-0.10364

0.037157

Модифицированный метод Эйлера:

X

ycv

zcv

y

z

y теор

z теор

y-y теор

0

1

0

1

0

1

0

0

0.2

1

-0.2

1

-0.18

0.983685

-0.14622

0.016315

0.4

0.96

-0.28

0.962

-0.244

0.947216

-0.20658

0.014784

0.6

0.904

-0.28

0.9096

-0.2314

0.905009

-0.20739

0.004591

0.8

0.848

-0.2288

0.85846

-0.17048

0.866913

-0.16826

0.008453

1

0.80224

-0.14688

0.818532

-0.08127

0.839397

-0.10364

0.020865

Схема Рунге – Кутта:

x

Y

z

k1

l1

k2

l2

k3

l3

k4

l4

0

1

0

0

-1

-0.1

-0.7

-0.07

-0.75

-0.15

-0.486

0.2

0.983667

-0.1462

-0.1462

-0.49127

-0.19533

-0.27839

-0.17404

-0.31606

-0.20941

-0.13004

0.4

0.947189

-0.20654

-0.20654

-0.13411

-0.21995

0.013367

-0.2052

-0.01479

-0.2095

0.112847

0.6

0.904977

-0.20734

-0.20734

0.10971

-0.19637

0.208502

-0.18649

0.187647

-0.16981

0.27195

0.8

0.866881

-0.16821

-0.16821

0.269542

-0.14126

0.332455

-0.13497

0.317177

-0.10478

0.369665

1

0.839366

-0.1036

-0.1036

0.367825

-0.06681

0.40462

-0.06313

0.393583

-0.02488

0.423019

Max(y-y теор)=4*10-5

Метод конечных разностей решения краевых задач для ОДУ

Постановка задачи: найти решение линейного дифференциального уравнения

, (1)

удовлетворяющего краевым условиям:. (2)

Теорема. Пусть . Тогда существует единственное решение поставленной задачи.

К данной задаче сводится, например, задача об определении прогибов балки, которая на концах опирается шарнирно.

Основные этапы метода конечных разностей:

1) область непрерывного изменения аргумента ([a,b]) заменяется дискретным множеством точек, называемых узлами: .

2) Искомая функция непрерывного аргумента x, приближенно заменяется функцией дискретного аргумента на заданной сетке, т.е. . Функция называется сеточной.

3) Исходное дифференциальное уравнение заменяется разностным уравнением относительно сеточной функции. Такая замена называется разностной аппроксимацией.

Таким образом, решение дифференциального уравнения сводится к отысканию значений сеточной функции в узлах сетки, которые находятся из решения алгебраических уравнений.

Аппроксимация производных.

Для аппроксимации (замены) первой производной можно воспользоваться формулами:

– правая разностная производная,

– левая разностная производная,

– центральная разностная производная.

т.е., возможно множество способов аппроксимации производной.

Все эти определения следуют из понятия производной как предела: .

Опираясь на разностную аппроксимацию первой производной можно построить разностную аппроксимацию второй производной:

(3)

Аналогично можно получить аппроксимации производных более высокого порядка.

Определение. Погрешностью аппроксимации n- ой производной называется разность: .

Для определения порядка аппроксимации используется разложение в ряд Тейлора.

Рассмотрим правую разностную аппроксимацию первой производной:

Т.е. правая разностная производная имеет первый по h порядок аппроксимации.

Аналогично и для левой разностной производной.

Центральная разностная производная имеет второй порядок аппроксимации.

Аппроксимация второй производной по формуле (3) также имеет второй порядок аппроксимации.

Для того чтобы аппроксимировать дифференциальное уравнение необходимо в нем заменить все производные их аппроксимациями. Рассмотрим задачу (1), (2) и заменим в(1) производные:

.

В результате получим:

(4)

Порядок аппроксимации исходной задачи равен 2, т.к. вторая и первая производные заменены с порядком 2, а остальные – точно.

Итак, вместо дифференциальных уравнений (1), (2) получена система линейных уравнений для определения в узлах сетки.

Схему можно представить в виде:

т.е., получили систему линейных уравнений с матрицей:

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

Решая полученную систему уравнений, мы получим решение исходной задачи.

Для решения таких СЛАУ имеется экономичный метод прогонки.

Рассмотрим метод прогонки для СЛАУ:

(1)

Решение данной системы ищем в виде:

(2)

Подставляя в первое уравнение, получим:

Здесь учтено, что данное соотношение должно выполняться при любом

Так как

, (3)

то подставляя (3) во второе уравнение, получим:

Сравнивая с (2) получим

.

Таким образом, можно найти все .

Тогда из последнего уравнения (1) находим:

Затем последовательно находим:

Таким образом, алгоритм метода прогонки можно представить в виде:

1) Находим

2) Для i=1,n-1: (4)

3) Находим

4) Для i=n-1 до 1 находим:

Шаги 1),2) – прямой ход метода прогонки, 3),4) – обратный ход метода прогонки.

Теорема. Пусть коэффициенты ai, bi системы уравнений при i =2, 3, …, n–1 отличны от нуля и пусть

при i =1, 2, 3, …, n. Тогда прогонка корректна и устойчива.

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

Для нашей краевой задачи имеем :

Тогда: , ,

Для нашей задачи условие устойчивости имеет вид:

.

Пусть

. (6)

Тогда

Пример. Найти решение задачи:

Выпишем разностную схему

Условие устойчивости примет вид

Возьмем .

Тогда

Или

Формулы прогонки были получены для СЛАУ (1):

Здесь x замены на u.

Следовательно,

Решим СЛАУ методом прогонки. Вычисления оформим в виде таблицы:

I

ai

ci

bi

fi

alfai

betai

ui

1

 

51

35

0.2

0.6863

-0.0039

0.4701

2

15

51

35

0.4

0.8598

-0.0113

0.6906

3

15

51

35

0.6

0.9186

-0.0202

0.8164

4

15

51

35

0.8

0.9403

-0.0296

0.9107

5

0

-1

 

1

 

 

1.0000

Порядок вычислений по формулам (4):

Ответ в столбце ui.

Если забыли формулы, то их можно легко вывести. Главное запомнить основную формулу:

Прямой ход

Обратный ход

На практике часто граничные условия могут иметь более общий вид.

Рассмотрим следующую краевую задачу:

Найти решение ОДУ 2-го порядка

,

удовлетворяющую краевым условиям:

В этом случае при построении разностной схемы необходимо еще аппроксимировать и краевые условия.

Аппроксимация:

В результате получим разностную схему:

Или

Мы получили СЛАУ типа (5) с трехдиагональной матрицей, решение которой также можно найти методом прогонки.


Метод Рунге Кутта, решение дифференциальных уравнений


    Для нахождения приближенного решения системы дифференциальных уравнений можно использовать метод Рунге — Кутта и метод Эйлера. [c.123]

    Для решения системы дифференциальных уравнений методом Рунге — Кутта необходимо выполнить вычисления в следующем порядке  [c.124]

    Численное интегрирование обыкновенны.х дифференциальных уравнений (задача Коши) выполняется одношаговыми методами, в которых решение в точке хп+ находится по известному решению в точке Хп- Наиболее распространенным одношаговым методом численного интегрирования является метод Рунге—Кутты четвертого порядка, и соответствии с которым решение уп л определяется по уп следующим образом  [c.147]

    Используют стандартную программу решения системы обыкновенных дифференциальных уравнений первого порядка (например, методом Рунге — Кутта или методом Адамса) с автоматическим выбором шага интегрирования в зависимости от требуемой точности вычисления. Эта программа позволяет определить значения концентрации х ( р, 0) и температуры совокупности точек, на которые разбивается интервал (О — Ь) интегрирования. [c.151]

    При численном решении систем дифференциальных уравнений наиболее часто используют методы Эйлера и Рунге — Кутта. С другими методами можно ознакомиться в книгах по вычислительной математике [2, 3]. Оба эти метода удобны при программировании решения на ЭВМ для тех случаев, когда все граничные (начальные) условия заданы при одном и том же значении аргумента. Охарактеризуем кратко эти методы. [c.145]

    Решение системы из двух дифференциальных уравнений методом Рунге-Кутта с фиксированным шагом [c.98]

    Нами выполнены расчеты результатов про-цесса по математическому описанию при тех же входных величинах, что и в промышленном аппарате с использованием стандартной программы решения системы дифференциальных уравнений методом Рунге—Кутта для ЭВМ М-20. Результаты расчетов при нескольких величинах ко показаны в табл. 8.2, где для удобства сравнения приведены и выходные опытные данные. При подборе ко в качестве исходного значения принята величина, рассчитанная на основании работы [147]. [c.181]

    Все кинетические константы, входящие в систему уравнений (6.5), были определены экспериментально. Решение системы дифференциальных уравнений осуществляли численным методом Рунге-Кутта 4-го порядка. [c.175]


    Оценим вреия решения задач математического моделирования целиком на ЦВМ с таким быстродействием. Ддя решения системы дифференциальных уравнений, например, методом Рунге-Кутта с заданной точностью обычно требуется [c.504]

    Решение полученной системы уравнений на ЭВМ при известных значениях Параметров К, Ту, Т , Та, /Со. о и заданных начальных условиях осуществляется по программе, в которую входит стандартная подпрограмма интегрирования дифференциальных уравнений по методам Рунге—Кутта, Хэмминга или Адамса [c.155]

    Задача расчета переходного процесса состоит в решении дифференциальных уравнений, описывающих состояние системы. При использовании цифровых вычислительных машин с этой целью применяют методы численного интегрирования дифференциальных уравнений. Достаточно широкое распространение при расчетах переходных процессов на ЭВМ получили методы Рунге— Кутта, Хэмминга и Адамса. Рассмотрим сущность этих методов на примере решения дифференциального уравнения первого порядка [c.154]

    Блок-схема алгоритма приведена в работе [36]. Для численного интегрирования системы обыкновенных дифференциальных уравнений, описывающих процесс каталитического риформинга, первоначально использовался метод Рунге—Кутта. Разработанная программа позволила эффективно интегрировать дифференциальные уравнения. Однако, как показала практика, на расчеты затрачивалось много времени. Для сокращения времени счета была составлена другая программа, использующая более быстрый метод Эйлера. Сравнение точности вычислений по этим двум методам решения системы дифференциальных уравнений приведено в таблице III. 2. Данные таблицы показывают, что [c.126]

    Решение системы линейных уравнений с учетом условий однозначности и соотношений, описывающих зависимость параметров процесса от искомых величин и, Т к х, возможно только численными методами, для чего дифференциальные уравнения в частных производных (6.101) — (6.103) записывались в конечно-разностном виде [33] по переменной координате /. Полученная система обыкновенных дифференциальных уравнений решалась методом Рунге — Кутта, для чего алгоритм расчета был реализован в виде ФОРТРАН-программы. [c.190]

    Поиск минимума функции Ф по переменным ка, Е осуществлялся градиентным методом. В процессе решения задачи было установлено, что Ф имеет овраги . Для движения по их дну применялся метод оврагов . Частные производные Ф(йо,, Ё) находились по разностной схеме. Решение системы дифференциальных уравнений (XI, 36), (XI. 37) производилось методом Рунге — Кутта с шагом, равным Vie объема реактора. Затраты машинного времени ЦВМ типа М-20 на поиск минимума составляли не менее 1,5—2 ч при достаточно хороших начальных приближениях. Минимальное значение Ф при использовании данных табл. XI. 4 и XI. 5 равно 4,2. [c.306]

    Решим дифференциальное уравнение (25) методом Рунге-Кутта с автоматическим выбором шага по программе,описанной в литературе, придавая константам уравнения различные значения из априорных соображений. Результаты решения указанного дифференциального уравнения приведены на рис. 7. [c.40]

    Эта математическая модель является упрощенной дополнительная подача пара в реактор-смеситель не учитывается, не рассматривается процесс кристаллизации гипса. Система нелинейных дифференциальных уравнений второго порядка (236), (257)—(259) аналитического решения не имеет. При реализации ма тематической модели реактора-смесителя поэтому были использованы численные методы интегрирования [359], в частности метод Рунге-Кутта четвертого порядка точности, дающий малую ошибку и легко программируемый (в нашем случае при моделировании на ЭВМ Наири-2 использовалась стандартная программа из библиотеки СП машины). [c.181]

    Возможны и более сложные реализации итерационных, рекуррентных и других подобных вычислений. Например, к ним сводится решение систем дифференциальных уравнений любыми разностными методами, например, Эйлера, Рунге—Кутта и др. [c.76]

    Сводка уравнений модели была дана в разд. 8.3. Для решения системы дифференциальных уравнений можно воспользоваться методами Рунге — Кутта или методами предсказания с коррекцией [55]  [c.211]

    Математическое описание модели [уравнения (10.5), (10.7) и (10.10)] приведено в разд 10.3. В гл. 8 уже упоминалось, что для решения системы обыкновенных дифференциальных уравнений могут быть использованы как метод Рунге — Кутта, так и методы предсказания с коррекцией. Из соображений, уже высказанных в гл. 8, был выбран первый из них. [c.240]

    Таким образом, система кинетических уравнений пиролиза углеводородов представляется в виде двух взаимосвязанных (через независимые переменные—компоненты реакционной смеси) систем — дифференциальных уравнений для молекулярных углеводородов (система Д) и алгебраических уравнений для радикалов (система А). Для решения каждой из этих систем могут быть использованы стандартные методы решения. В частности, для решения системы Д был использован метод Рунге — Кутта четвертого порядка с автоматическим выбором шага интегрирования, а для системы А — метод простой итерации [108, 109]. Решение системы А выполняется в тех же точках, что и для системы Д. [c.46]


    Программу расчета составим в соответствии с укрупненной структурой ал горитма, показанной на рис. 5.17. При этом воспользуемся стандартной под программой ЯКОЗ решения дифференциальных уравнений по методу Рунге— Кутта. В подпрограмме предусмотрено вычисление весовых коэффици12нтов значения которых могут быть взяты одинаковыми при одинаковом порядке оП ределяемых при расчете переменных,. Для выполнения этого условия перейдем [c.387]

    Биотехнологические процессы моделируются с помощью систем дифференциальных уравнений. Только немногие из этих систем можно решить аналитически, чаще требуется применение численных методов интегрирования прикидочное решение достигается с помощью метода Эйлера, более точное решение дает метод Рунге Кутта. [c.57]

    Программа моделирования на цифровой ЭВМ. Программу моделирования реактора на цифровой ЭВМ применяли для интегрирования уравнений материального и теплового баланса реактора идеального вытеонения. Численные решения системы нелинейных дифференциальных уравнений получали методом Рунге-Кутта четвертого порядка. Всю систему дифференциальных уравнений интегрировали по длине реактора и получали концентрационные и температурные профили. Основная программа была управляющей, а уравнения скорости реакций и термодинамические характеристики вычисляли в подпрограмме 5иЬги11пе. В этой подпрограмме реализуется печать результатов каждого шага интегрирования, содержащих информацию по составу и температуре. Кроме того, рассчитывали и печатали значения выходов, селективностей и степеней превращения. Таким образом, имелась подробная информация по ходу моделирования для широких диапазонов изученных условий. [c.292]

    Алгебраические уравнения (I) – (2) и описание слоя идеального вытеснения (3) – (4) объединим в первую группу. Путем сведения к нестационарной задаче трансцендентные уравнения (I) -(2), непосредственное решение которых встречает определенные трудности, преобразуются в систему обыкновенных дифференциальных уравнений, как (3) – (4). Эти уравнения представляют собой задачу Коши, для решения которой имеется ряд хорошо разработанных методов, например, Рунге-Кутта, Симпсона и др. Отметим, что для (I) – (2) представляет интерес решение стационарной задачи, т.е. при (если t – переменная, играющая роль времени [c.136]

    Существенным моментом при создании специализированных пакетов прикладшхх программ является использование одного или ограниченной совокупности методов для решения широкого класса задач. Значительный опыт по разработке таких систем накоплен при решении дифференциальных уравнений, для описания динамических систем (расчет траекторий полета спутников, баллистика и т. д.). К таким системам можно отнести системы MIDAS [17], MIMI [18], в основе которых используются формулы Рунге— Кутта различного порядка. [c.275]

    Результаты расчетов. Численное решение полученной системы уравнений осуществляется на основе комбинации явного (метод Рунге — Кутта) и нолунеявного (метод Михельсона) методов решения обыкновенных дифференциальных уравнений. Размерность системы определяется дифференциальными уравнениями, описывающими как непрерывную (17)—(20), так и дискретную (21) —(23) фазы для каждого класса капель. По мере исчезновения г-го класса размерность уменьшается на число уравнений, описывающих его. [c.77]

    Матрица, содержащая таблицу значений решения задачи Кощи на интервале от х 1 до х2 для системы обыкновенных дифференциальных уравнений, вычисленную методом Рунге-Кутта с переменным щагом и начальными условиями в векторе v, причем правые части системы записаны в D, п — число шагов, к — максимальное число промежуточных точек решения, s — минимально допустимый интервал между точками (только для Malh ad Professional) [c.452]

    Матрица решения системы обыкновенных дифференциальных уравнений численным методом Рунге—Ку гга на интервале от х1 до х2 с переменным шагом, при минимальном числе шагов п, причем правые части уравнений в символьной форме задаются в векторе D, а начальные условия — в векторе V (только для Math ad Professional) Матрица решения системы обыкновенных дифференциальных уравнений методом Рунге—Кутта на интервале от х1 до х2 при фиксированном числе шагов п, причем правые части уравнений записаны в символьном векторе D, а начальные условия — в векторе v [c.453]

    Моделирование процесса эмульсионной полимеризации на ЦВМ. Для численного решения задачи (3.47)—(3.63) с начальными, граничными условиями и условиями сопряжения (3.64) — (3.68) система дифференциальных уравнений приводилась к безразмерному виду и решалась методом прямых с применениеи процедуры Рунге—Кутта—Мерсона на ЦВМ Минск-32 . [c.156]

    В ранних работах, выполненных методом классических траекторий, очень популярным был метод Рунге—Кутта—Гилла 4-го порядка [265]. Позже стали применяться многошаговые методы высокого (до 16-го) порядка точности, в основном использовалась процедура Адамса-Мултона 4-го порядка [92], метод экстраполяций [219]. В большинстве работ численное решение систем дифференциальных уравнений осуществлялось с постоянным шагом интегрирования. В [66] было проведено сравнение эффективности различных процедур численного интегрирования для решения систем обыкновенных дифференциальных уравнений. [c.77]

    Чтобы количественный анализ динамики объемного привода был достоверным, необходимо учитывать зависимость коэффициентов уравнений (2.159) от переменных величин. В этом случае дифференциальные уравнения будут нелинейными и для их -решения используют численные методы интегрироьв-ния [17]. Эти методы трудоемки и требуют применения ЭВМ. Разработаны программы численного интегрирования обыкновенных дифференциальных уравнений на основе неяоных методов Рунге-Кутта, Адамса и др., которые обеспечивают автоматический выбор временного шага и порядка метода. [c.150]

    Программа МАТ21 дает решение обыкновенного дифференциального уравнения первой степени методом Рунге — Кутта. [c.243]

    У1 ( А + 1) У ( а) б г 4, г ( 2, г + 3, /)) = П. Программа М.АТ22 дает решение системы обыкновенных дифференциальных уравнений первой степени методом Рунге — Кутта. [c.245]

    Рассмотрим подробно решение обратной задачи в среде MS Ex el. Для решения задачи Коши системы двух дифференциальных уравнений используем метод Рунге—Кутта четвертого порядка. Расчеты будем проводить по формулам  [c.275]

    Если попытаться решить жесткую систему дифференш1альных уравнений, описываюшую вполне реальную химическую реакцию, обычными численными методами, например широко распространенным в естественных науках методом Рунге — Кутта, то, как правило, получаются ложные результаты. Если интегрировать дифференциальные уравнения с очень малым шагом, то даже при больших затратах машинного времени решение соответствует незначительным изменениям в системе, т. е. очень малой степени преврашения. [c.395]

    Предположим, что х t) и х t) заданы ординатами х (ti), х ti) при i = 1, 2,. . d. Тогда для нахождения Ф по формуле (V-30) требуется d раз вычислять значения функции /. Если же в Ф входит решение дифференциального уравнения, найденное, например, методом Рунге — Кутта порядка q (обычно q — А), то для определения значения Ф (а) требуется дТIh раз вычислять функцию / (здесь h — шаг численного интегрирования). Обычно значение Tlh не менее чем на порядок больше числа d, поэтому затраты машинного времени на вычисление Ф по формуле (V-30) сокращаются примерно в 40 раз. [c.295]

    Система уравнений решалась на цифровой вычислительной маш не при использовании метода Рунге-Кутта. Для решения дифференциальных уравнений величина Ко( е)п, определяющая концентрацию в органической фазе в стационарном состоянии, определялась экспериментально. Причем, (/Сса )п приближалась к Коар) эксп только с увеличением числа ячеек. [c.138]

    Задаваясь каким-либо значением а= у с учетом ограничения (V.12), по формуле (V.8) вычисляем управление U x), соответствующее av. Подставляем полученное U x) в систему дифференциальных уравнений (V.1) и решаем полученную систему при начальных условиях (V.2). Решение осуществим на ЭВМ с помощью метода Рунге — Кутта при автоматическом выборе длины шага интегрирования. В результате решения найдем Xi(L), Х Ь), Xs(L), X4(L) из формулы (V.3) следует, 4To7 i/ = / (aJ. Таким [c.108]

    Седьмая глава является одной из основных глав книги. Здесь на примерах теплообменника и ректификационной колонны показана методика использования численных методов решения задач. Авторы связывают расчетные параметры с изучаемым процессом. Рассматриваются методы преобразования дифференциальных уравнений в частных производных в систему обыкновенных дифференциальных уравнений и затем в разностные уравнения. Сравнивается использование различных методов (Эйлера, Рунге—Кутта, Крэнка—Никольсона и метода авторов) с точки зрения сходимости, точности и возможности расчета с помощью цифровых вычислительных машин (ЦВМ). Приводится расчет многокомпонентной ректификационной колонны. В заключение дается обзор численных методов. Следует отметить, что опущены некоторые математические рассуждения, очень простые для математиков и необходимые для понимания химикам-технологам. [c.7]

    В математическом отношении расчет периодической ректификации многокомпопентной смеси в приближении теоретической тарелки сводится к интегрированию обширной системы обык]к )вениых дифференциальных уравнений. На практике, главным образом, используются два метода численого решения задачи Коши машинные варианты метода Рунге—Кутта [1, 2] и неявный одношаговый конечно-разностный метод, имеющий в основе квадратурную формулу трапеций [3, 4]. В первом случае известные трудности представляет нахождение явного вида прои родной от температуры по времени, кроме того, система уравнений периодической ректификации относится к типу жестки.х систем, для которых методы Рунге—Кутта могут потребовать очень малого шага интегрирования или вообще ие будут работать [5]. Неявный метод более подходит для интегрирования жест.ких систем, но требуег большего объема вычислений иа каждом шаге, поскольку сводит решение нестационарной задачи к последовательному решению нелинейных систем алгебраических уравнений. [c.62]

    Построена упрощенная математическая модель кинетики фотохимического синтеза четырех основных изомеров гексахлорциклогексана в изотермических условиях. Модель представляет собой систему из 20 лилейных дифференциальных уравнений псевдопервого порядка. Решение системы уравнений проводи лось численным методом Рунге—Кутта на ЭВМ. Полученная модель удовлетворительно аппроксилирует опытные да1тые в области, представляющей практический интерес. [c.123]

    I. Вычисление характеристик периодического процесса по уравнениям типа (1.36), (1.37) и т.п. (в зависимости от вида кинетического модуля) с использованием методов численного решения систем обыкновенных дифференциальных уравнений (Рунге — Кутта, Адамса или других из имеющихся в системе математического обеспечения используемой ЭВМ). Отметим, что в качестве начальных условий при решении по уравнениям (1У.2) должны быть выбраны значения выхода предыдущего реактора [c.140]

    О численном решении задачи об автоволне в нестационарной постановке. В качестве метода решения начально-краевой задачи (1.8) был выбран метод прямых [29]. Использовалась аппроксимация второго порядка точности по пространственной переменной. Полученная система обыкновенных дифференциальных уравнений интегрировалась с помош,ью неявной схемы Рунге-Кутта 5-го порядка точности, программная реализация которой эффективно учитывала ленточный вид матрицы Якоби правых частей. [c.64]

    Заканчивая данный раздел, сделаем некоторые замечания относительно использования явных методов для численного решения жестких систем дифференциальных уравнений. В ряде случаев возникает необходимость применения явных формул для решения жестких задач. Это требуется, например, при большой размерности дифференциальной задачи. Алгоритмы на основе неявных или полуявных формул, как правило, используют обращение матрицы Якоби, что в данном случае есть отдельная трудновыполнимая задача. В такой ситуации предпочтительнее использовать алгоритмы на основе явных формул, если жесткость задачи позволяет за разумное время получить приближение к решению. Современные алгоритмы на основе явных формул в большинстве своем не приспособлены для решения жестких задач по следующей причине. Обычно алгоритм управления величиной шага строится на контроле точности численной схемы. Это естественно, так как основным критерием является точность вычисления решения. Однако при применении алгоритмов интегрирования на основе явных формул для решения жестких задач этот подход приводит к потере эффективности и надежности, ибо вследствие противоречивости требований точности и устойчивости шаг интегрирования раскачивается, что приводит либо к большому количеству возвратов (повторных вычислений решения), либо к АВОСТу. Этого можно избежать, если наряду с точностью контролировать и устойчивость численной схемы. Может быть предложен способ контроля устойчивости явных методов и алгоритм интегрирования с контролем точности и устойчивости на основе явной формулы типа Рунге—Кутта второго порядка точности  [c.279]


численные методы решения ОДУ / Хабр

Очень кратко рассмотрев основы механики в

предыдущей статье

, перейдем к практике, ибо даже той краткой теории что была рассмотрена хватит с головой.

Итак, задача:

Камень бросают вертикально, без начальной скорости с высоты h = 100 м. Пренебрегая сопротивлением воздуха, определить закон движения камня, как функцию высоты камня над поверхностью Земли от времени. Ускорение свободного падения принять равным 10 м/с2

Простая задачка? Да элементарная, имеющее аналитическое решение, которое легко напишет мало-мальски грамотный школьник. Но эта простая задача послужит нам весьма показательным примером



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

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

Следующим допущением будет полагать что Земля плоская, так как высота, с которой мы бросаем камень пренебрежимо мала в сравнении с размерами планеты, а значит кривизной её поверхности можно пренебречь. Тогда сила тяжести может считаться постоянной, направленной перпендикулярно поверхности Земли и равной по модулю

где g — ускорение свободного падения у поверхности Земли. Теперь пришло время составить уравнения движения камня. Помните эти уравнения

Левая их часть нас пока не интересует, а вот в правой стоят суммы проекций сил, приложенных к точке на оси координат. Пусть оси x и y располагаются на поверхности, а ось z направлена вверх перпендикулярно ей. Сила одна единственная, её проекции на оси x и y равны нулю, а на ось z проекция отрицательна, так как сила направлена против направления оси, то есть

Масса камня, очевидно не равна нулю, значит можно спокойно поделить обе части получившихся уравнений на неё

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

То что у нас получилось не много не мало — математическая модель процесса происходящего в задаче. Пафосно, да?

Нет. Анализируя это уравнение мы делаем вывод, например, что масса камня не оказывает влияния на закон его движения, ведь массы в этом уравнении нет. Видите, даже не решив уравнения, мы уже формально доказали справедливость опыта с пером и кусочком свинца в вакууме, который любят показывать в школе (а некоторые повторили его на Луне).

Аналитическое решение получить просто, даже не буду заморачиваться, оно такое

А вот как решить это численно? И что это вообще такое — «численно»?

Какой такой первый порядок? Я же говорил в прошлый раз, что уравнения движения имеют второй порядок. Всё правильно, но большинство методов решения диффур на компьютере умеют решать только уравнения первого порядка. Есть методы прямого интегрирования уравнений второго порядка (например метод Верле), но о них не сейчас.

Во-первых, это уравнение относится к такому типу, что допускает понижение порядка. Правая часть не зависит от неизвестной функции (там нет z), поэтому вспоминаем, что

проекция ускорения на ось z равна первой производной проекции скорости на ту же оcь z. Ну классно, тогда

вот вам и уравнение первого порядка. Не всегда этот номер проходит (не буду я сейчас про форму Коши!), но в данном случае всё в порядке. Будем искать не координату а скорость точки. Что дальше-то? А дальше

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

Что получается? А вот что

Мы получили приращение скорости. Отрицательное приращение. Как это так, камень, падая вниз будет разгонятся же! Да, будет. Его скорость, вектор его скорости, будет направлен вниз. А значит проекция этого вектора на ось z будет отрицательной. Всё правильно, мы получаем растущую по абсолютному значению проекцию вектора, направленного вниз. Мы знаем, начальное значение скорости — ноль, а значит

Пользуясь тем что мы можем вычислить приращение скорости, посчитаем, какова будет скорость скажем через 0.1 секунды

а ещё через 0.1 секунды

и ещё через 0.1 секунды

Хм, так мы можем продолжать довольно долго, но ограничимся промежутком времени в одну секунду

Время, с Скорость, м/с
0.0 0.0
0.1 -1.0
0.2 -2.0
0.3 -3.0
0.4 -4.0
0.5 -5.0
0.6 -6.0
0.7 -7.0
0.8 -8.0
0.9 -9.0
1.0 -10.0

То есть, воспользовавшись формулой

мы получили зависимость скорости точки от времени. А всего-то нужно взять значение скорости в текущий момент времени, и добавить к нему приращение, которое скорость получит в новый момент времени, отстоящий от текущего на секунд. Приращение времени называется здесь шагом интегрирования. А приращение мы вычисляем как значение производной искомой функции в текущий момент времени умноженное на шаг. Просто? Да просто конечно. И та формула, которую я написал, имеет название название — явный метод Эйлера для численного решения дифференциальных уравнений. Это так называемая рекуррентная формула, когда новое значение вычисляемой величины зависит от её предыдущего значения.

А что же с высотой точки над Землей? Да всё аналогично, смотрите

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

и по этой формуле добавим в нашу таблицу ещё одну колонку

и так далее

Время, с Скорость, м/с Высота, м
0.0 0.0 100.0
0.1 -1.0 100.0
0.2 -2.0 99.9
0.3 -3.0 99.7
0.4 -4.0 99.4
0.5 -5.0 99.0
0.6 -6.0 98.5
0.7 -7.0 97.9
0.8 -8.0 97.2
0.9 -9.0 96.4
1.0 -10.0 95.5

Хм, ну, во-первых, заметно, что высота меняется у нас уже неравномерно, так как скорость со временем меняется. Теперь наша производная сама зависит от времени. Но уже на первом шаге, мы замечаем неладное — скорость уже есть, а вот высота по прежнему 100 метров. Как так?

Это вышло потому, что на каждом шаге мы полагаем производную (скорость) постоянной. Метод не дает информации о том, что происходит с решением между шагами. Соответственно накапливается погрешность, сравним полученное решение с точным

Время, с Скорость, м/с Высота, м Точное решение, м
0.0 0.0 100.0 100.00
0.1 -1.0 100.0 99.95
0.2 -2.0 99.9 99.80
0.3 -3.0 99.7 99.55
0.4 -4.0 99.4 99.20
0.5 -5.0 99.0 98.75
0.6 -6.0 98.5 98.20
0.7 -7.0 97.9 97.55
0.8 -8.0 97.2 96.80
0.9 -9.0 96.4 95.95
1.0 -10.0 95.5 95.00

Да, наш камень как будто зависает в воздухе. Численное решение отстает от аналитического, и чем дальше мы считаем, тем выше погрешность счета. Погрешность накапливается, так как на каждом шаге мы берем всё более и более грубое приближение. Что делать?

Во первых, можно уменьшить шаг. Скажем в 10 раз, пусть секунды

Время, с Скорость, м/с Высота, м Точное решение, м
0.0 0.0 100.0 100.00
0.1 -1.0 99.96 99.95
0.2 -2.0 99.81 99.80
0.3 -3.0 99.57 99.55
0.4 -4.0 99.22 99.20
0.5 -5.0 98.78 98.75
0.6 -6.0 98.23 98.20
0.7 -7.0 97.59 97.55
0.8 -8.0 96.84 96.80
0.9 -9.0 96.00 95.95
1.0 -10.0 95.05 95.00

Уже лучше, погрешность в конце счета не превышает 0,05 метров, и это в 10 раз меньше предыдущего значения. Можно предположить, что уменьшив шаг ещё в 10 раз мы получим ещё более точное решение. Я схитрил, выводя значения только для 10 точек с шагом 0.1, на самом деле, чтобы получить такую таблицу нужны уже 100 итерации а не 10. При шаге 0.001 потребуется уже тысяча итераций, а результат будет таким

Время, с Скорость, м/с Высота, м Точное решение, м
0.0 0.0 100.0 100.00
0.1 -1.0 99.9505 99.95
0.2 -2.0 99.8010 99.80
0.3 -3.0 99.5515 99.55
0.4 -4.0 99.2020 99.20
0.5 -5.0 98.7525 98.75
0.6 -6.0 98.2030 98.20
0.7 -7.0 97.5535 97.55
0.8 -8.0 96.8040 96.80
0.9 -9.0 95.9545 95.95
1.0 -10.0 95.0050 95.00

Если вы попробовали выполнить эти расчеты в ручную, то понимаете теперь насколько они однообразны и трудоемки, если нужна высокая точность. Именно поэтому расцвет численного моделирования совпал с появлением компьютеров. Они как раз и нужны для того, чтобы быстро выполнять множество однообразных операций над числами.

Метод Эйлера самый простой из известных методов интегрирования дифференциальных уравнений. Из нашего простого примера видно, что погрешность метода прямо пропорциональна шагу интегрирования, и это действительно так. Такие методы называются методами 1-го порядка точности.

Точность расчетов даже на шаге 0.1 можно улучшить, если мы применим, скажем метод Рунге-Кутты 4-го порядка точности. Но это отдельная история.

Задумайтесь… Мы рассмотрели очень простой пример. Мы даже не применяли компьютер, но уже понимаем принцип, по которому работают те самые мощные в мире суперкомпьютеры, что моделируют ранние этапы жизни Вселенной. Конечно, там всё устроено гораздо сложнее, но принцип лежит этот же самый.

Представьте себе, какой мощный инструмент вы получаете в свои руки. Эта последняя статья, где мы не будем применять компьютер. Я обещал Octave. В следующий раз будет именно он.

Редукция порядка

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

Тип 1: уравнения второго порядка с отсутствующей зависимой переменной

Тип 2: нелинейные уравнения второго порядка с отсутствующей независимой переменной

Тип 3: однородные линейные уравнения второго порядка, в которых известно одно (ненулевое) решение

Тип 1: уравнения второго порядка с отсутствующей зависимой переменной. Примеры таких уравнений включают

Определяющая характеристика заключается в следующем: зависимая переменная y явно не появляется в уравнении. Этот тип уравнения второго порядка легко сводится к уравнению первого порядка преобразованием

Эта замена, очевидно, подразумевает y ″ = w ′, и исходное уравнение становится уравнением первого порядка для w . Решить для функции w ; затем интегрируйте его, чтобы получить y .

Пример 1 : Решите дифференциальное уравнение y ′ + y ″ = w .

Поскольку зависимая переменная y отсутствует, пусть y ′ = w и y ″ = w ′. Эти замены преобразуют данное уравнение второго порядка в уравнение первого порядка

в стандартной форме. Применяя метод решения таких уравнений, сначала определяется интегрирующий коэффициент,

, а затем умножили обе части уравнения, получив

Следовательно,

Теперь, чтобы получить решение y исходного уравнения второго порядка, проинтегрируйте:

Это дает

Ссылаясь на теорему B, обратите внимание, что это решение подразумевает, что y = c 1 e x + c 2 является общим решением соответствующего однородного уравнения и что y = ½ x 2 x является частным решением неоднородного уравнения.(Это конкретное дифференциальное уравнение можно было бы также решить, применив метод решения линейных уравнений второго порядка с постоянными коэффициентами.)

Пример 2 : Решите дифференциальное уравнение

Опять же, зависимая переменная y отсутствует в этом уравнении второго порядка, поэтому ее порядок будет уменьшен путем замены y ′ = w и y ″ = w ′:

, который можно записать в стандартной форме

Интегрирующий коэффициент здесь

.

, которое используется для умножения обеих частей уравнения, что дает

Интегрирование w дает y :

Положив c 1 = ⅓ c 1 , общее решение можно записать

Пример 3 : Набросок решения IVP

Хотя это уравнение является нелинейным [из-за члена ( y ′) 2 ; ни y , ни любая из его производных не могут быть возведены в какую-либо степень (кроме 1) в линейном уравнении], замены y ′ = w и y ″ = w ′ все равно не будут свести это к уравнению первого порядка, поскольку переменная y не появляется явно.Дифференциальное уравнение преобразуется в

который отделяется:

Поскольку y ′ = w , интегрирование дает

Теперь примените начальные условия для определения констант c 1 и c 2 .

Поскольку c 1 = 1, из первого условия также следует c 2 = 1.Таким образом, решение этой IVP (по крайней мере, для x > −1) составляет

, график которого показан на рис.


Рисунок 1

Тип 2: нелинейные уравнения второго порядка с отсутствующей независимой переменной. Вот пример такого уравнения:

Определяющая характеристика такова: независимая переменная x явно не фигурирует в уравнении.

Метод уменьшения порядка этих уравнений второго порядка начинается с той же замены, что и для уравнений типа 1, а именно с замены y ′ на w . Но вместо того, чтобы просто записать y ″ как w ′, уловка здесь состоит в том, чтобы выразить y ″ через первую производную по y . Это достигается с помощью цепного правила:

Следовательно,

Эта замена вместе с y ′ = w сведет уравнение типа 2 к уравнению первого порядка для w .Как только w определено, проинтегрируйте, чтобы найти y .

Пример 4 : Решите дифференциальное уравнение

Замены y ′ = w и y ″ = w ( dw / dy ) преобразуют это уравнение второго порядка для y в следующее уравнение первого порядка для w :

Следовательно,

Выражение w = 0 означает y ′ = 0, и, таким образом, y = c является решением для любой константы c .Второе утверждение представляет собой разделимое уравнение, и его решение происходит следующим образом:

Теперь, поскольку w = dy / dx , этот последний результат становится

, который дает общее решение, неявно выраженное следующим образом:

Следовательно, полное решение данного дифференциального уравнения равно

Тип 3: однородные линейные уравнения второго порядка, для которых известно одно (отличное от нуля) решение. Иногда можно определить решение дифференциального уравнения второго порядка путем проверки, что обычно сводится к успешным пробам и ошибкам с помощью нескольких особенно простых функций. Например, вы можете обнаружить, что простая функция y = x является решением уравнения

или что y = e x удовлетворяет уравнению

Конечно, метод проб и ошибок – не лучший способ решить уравнение, но если вам повезет (или вы достаточно натренированы), чтобы действительно найти решение путем проверки, вы должны быть вознаграждены.

Если известно одно (ненулевое) решение однородного уравнения второго порядка, существует простой процесс определения второго, линейно независимого решения, которое затем может быть объединено с первым для получения общего решения. Пусть y 1 обозначает функцию, которая, как вы знаете, является решением. Тогда пусть y = y 1 v ( x ), где v – функция (пока неизвестная). Подставляем y = y 1 v в дифференциальное уравнение и получаем уравнение второго порядка для v .Это окажется уравнение типа 1 для v (поскольку зависимая переменная v не будет отображаться явно). Используйте метод, описанный ранее, чтобы найти функцию v ; затем подставьте в выражение y = y 1 v , чтобы получить желаемое второе решение.

Пример 5 : Дайте общее решение дифференциального уравнения

Как упоминалось выше, легко найти простое решение y = x .Обозначив это известное решение как y 1 , подставим y = y 1 v = xv в данное дифференциальное уравнение и решим относительно v . Если y = xv , то производные равны

Подстановка в дифференциальное уравнение дает

Обратите внимание, что это результирующее уравнение является уравнением типа 1 для v (поскольку зависимая переменная v не отображается явно).Итак, если предположить, что v ′ = w и v ”= w ′, это уравнение второго порядка для v становится следующим уравнением первого порядка для w :

Интегрирующий коэффициент для этого стандартного линейного уравнения первого порядка равен

.

и умножение обеих частей (*) на μ = x дает

Игнорировать константу c и интегрировать для восстановления v :

Умножьте это на y 1 , чтобы получить желаемое второе решение,

Общее решение исходного уравнения представляет собой любую линейную комбинацию y 1 = x и y 2 = x In | x |:

Это согласуется с общим решением, которое было бы найдено, если бы эта проблема была решена с использованием метода решения равноразмерного уравнения.

Пример 6 : Определите общее решение следующего дифференциального уравнения, учитывая, что ему удовлетворяет функция y = e x :

Обозначая известное решение как y 1 , подставляем y = y 1 v ′ = e x v в дифференциальное уравнение. Если y = e x u , производные равны

.

Подстановка в данное дифференциальное уравнение дает

, которое упрощается до следующего уравнения второго порядка типа 1 для v :

Положив v ′ = w , затем переписав уравнение в стандартной форме, получим

Интегрирующий коэффициент в данном случае равен

.

Умножение обеих частей (*) на μ μ = e x / x дает

Игнорировать константу c и интегрировать для восстановления v :

Умножьте это на y 1 , чтобы получить желаемое второе решение,

Общее решение исходного уравнения представляет собой любую линейную комбинацию y 1 и y 2

Произошла ошибка при настройке вашего пользовательского файла cookie

Произошла ошибка при настройке вашего пользовательского файла cookie

Этот сайт использует файлы cookie для повышения производительности.Если ваш браузер не принимает файлы cookie, вы не можете просматривать этот сайт.

Настройка вашего браузера для приема файлов cookie

Существует множество причин, по которым cookie не может быть установлен правильно. Ниже приведены наиболее частые причины:

  • В вашем браузере отключены файлы cookie. Вам необходимо сбросить настройки своего браузера, чтобы он принимал файлы cookie, или чтобы спросить вас, хотите ли вы принимать файлы cookie.
  • Ваш браузер спрашивает вас, хотите ли вы принимать файлы cookie, и вы отказались.Чтобы принять файлы cookie с этого сайта, используйте кнопку «Назад» и примите файлы cookie.
  • Ваш браузер не поддерживает файлы cookie. Если вы подозреваете это, попробуйте другой браузер.
  • Дата на вашем компьютере в прошлом. Если часы вашего компьютера показывают дату до 1 января 1970 г., браузер автоматически забудет файл cookie. Чтобы исправить это, установите правильное время и дату на своем компьютере.
  • Вы установили приложение, которое отслеживает или блокирует установку файлов cookie.Вы должны отключить приложение при входе в систему или проконсультироваться с системным администратором.

Почему этому сайту требуются файлы cookie?

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

Что сохраняется в файле cookie?

Этот сайт не хранит ничего, кроме автоматически сгенерированного идентификатора сеанса в cookie; никакая другая информация не фиксируется.

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

Дифференциальные уравнения – метод Эйлера – малый шаг

Рассмотрим линейное дифференциальное уравнение следующего вида: y ′ = dydx = f (x, y). y ‘= \ dfrac {dy} {dx} = f (x, y) .y ′ = dxdy = f (x, y). При решении дифференциального уравнения мы обычно сталкиваемся с уравнением, которое можно решить с помощью определенных методов, но в большинстве случаев дифференциальные уравнения не могут быть представлены в упрощенной форме. Чтобы обойти этот барьер, был разработан численный метод приближения. Один из самых простых и старых методов аппроксимации дифференциальных уравнений. известен как метод Эйлера .Метод Эйлера – это метод первого порядка, что означает, что локальная ошибка пропорциональна квадрату размера шага, а глобальная ошибка пропорциональна размеру шага. Метод Эйлера часто служит основой для построения более сложных методов.

Метод Эйлера основан на том факте, что рядом с точкой функция и ее касательная имеют почти одинаковое значение. Пусть hhh будет инкрементным изменением координаты xxx, также известным как размер шага.

ди-джеи

Из рисунка выше у нас есть наклон касательной в точке (x0, y0) (x_ {0}, y_ {0}) (x0, y0) как y ′ (x0) = f (x0 , y0) y ‘(x_ {0}) = f (x_ {0}, y_ {0}) y ′ (x0) = f (x0, y0).Учитывая наклон касательной и начальную точку (x0, y0) (x_ {0}, y_ {0}) (x0, y0), мы хотим узнать значение y1y_ {1} y1, расположенное при x1 = x0 + hx_ {1} = x_ {0} + hx1 = x0 + h.

Имеем из координатной геометрии

y1 = y0 + hf (x0, y0) .y_ {1} = y_ {0} + hf (x_ {0}, y_ {0}). Y1 = y0 + hf (x0, y0).

Теперь мы рекурсивно продолжаем с (x1, y1) (x_ {1}, y_ {1}) (x1, y1), чтобы найти значения (x2, y2) (x_ {2}, y_ {2}) (х2, у2). Наклон касательной в точке (x1, y1) (x_ {1}, y_ {1}) (x1, y1) равен y ′ (x1) = f (x1, y1) y ‘(x_ { 1}) = f (x_ {1}, y_ {1}) y ′ (x1) = f (x1, y1).Затем, используя тот же шаг, что и выше, у нас есть

y2 = y1 + hf (x1, y1) .y_ {2} = y_ {1} + hf (x_ {1}, y_ {1}). Y2 = y1 + hf (x1, y1).

В общем виде имеем

yn = yn − 1 + hf (xn − 1, yn − 1) .y_ {n} = y_ {n-1} + hf (x_ {n-1}, y_ {n-1}). Yn = yn − 1 + hf (xn − 1, yn − 1).

Теперь, чтобы понять, для чего мы только что сделали, давайте посмотрим на рисунок ниже, который показывает, что происходит, когда ччч становится все меньше и меньше:

fj

Синяя линия – это линия с наименьшим значением hhh, поэтому по мере того, как hhh становится все меньше и меньше, значение тангенса становится все ближе и ближе к реальному значению.{2} + 3xf ′ (x) = x2 + 3x. Используя метод Эйлера с шагом 1,1,1, найти полученное приближение f (5) .f (5) .f (5).


Используя метод Эйлера с шагом 1,1,1, имеем

f (3) = f (2) + 1 × f ′ (2) = 20f (4) = f (3) + 1 × f ′ (3) = 20 + 18 = 38f (5) = f (4) + 1 × f ′ (4) = 38 + 40 = 78. {2} +2 (1)] (0.{2} +2 (7,5)] (0,5) \\ & = 990,75. \ _\квадрат \ end {align} f (1.5) f (2) f (2.5) f (8) = f (1) + f ′ (1) × 0,5 = 2 + [6 (1) 2 + 2 (1)] (0,5) = 2 + 4 = 6 = f (1,5) + f ′ (1,5) × 0,5 = 6 + [6 (1,5) 2 + 2 (1,5)] (0,5) = 14,25 = f (2) + f ′ (2) × 0,5 = 14,25 + [6 (2) 2 + 2 (2)] (0,5) = 28,25⋅⋅ = f (7,5) + f ′ (7,5) × 0,5 = 814,5 + [6 (7,5) 2+ 2 (7,5)] (0,5) = 990,75. □

Экстраполяция методов первого порядка для параболических дифференциальных уравнений с частными производными, II в JSTOR

Абстрактный

В Части I этой статьи (SIAM J. Numer. Anal., 15 (1978), стр.1212-1224), новые алгоритмы были введены для решения параболических дифференциальных уравнений, в которых высокочастотные компоненты присутствовали в решении и для которых A 0 -устойчивые методы, представленные классическим методом Кранка-Николсона, были менее чем удовлетворительными. Представленные алгоритмы были основаны на простой экстраполяции обратного метода Эйлера, которая дала L 0 -стабильность. Во всех случаях представленные алгоритмы имели второй порядок точности по времени. В данной статье некоторые из предыдущих алгоритмов обобщены до более высоких порядков точности.В частности, дается более подробное рассмотрение алгоритма второго порядка, который естественным образом приводит к методам третьего и четвертого порядков. Новые алгоритмы тестируются на уравнении теплопроводности с постоянными коэффициентами, в котором существует разрыв между начальными и граничными значениями.

Информация о журнале

SIAM Journal on Numerical Analysis содержит исследовательские статьи по разработке и анализу численных методов, в том числе их сходимость, стабильность и анализ ошибок, а также связанные результаты в функциональный анализ и теория приближений.Вычислительные эксперименты также включены новые типы числовых приложений.

Информация об издателе

“Общество промышленной и прикладной математики является ведущим международная ассоциация прикладной математики и ее публикации мог бы стать ядром адекватного собрания по математике. Один из Целями этой организации является обеспечение обмена информацией между университет и промышленность стали более гладкими. Он превосходно выполняет эту задачу и многие из ведущих академических институтов мира являются членами.” – Журналы для библиотек, восьмое издание, 1995, Р. Р. Боукер, Нью-Провиденс, Нью-Джерси Общество промышленной и прикладной математики (SIAM), штаб-квартира в Филадельфии, была основана в 1951 году для продвижения применения математики в науку и промышленность, продвигать математические исследования и предоставлять средства массовой информации для обмена информацией и идеями между математики, инженеры и ученые. SIAM имеет обширную программу публикаций в прикладных и вычислительных математика, в том числе 11 престижных исследовательских журналов.Для полного описание наших журналов и недавно анонсированных SIAM Journals Online, доступ http://www.siam.org/.

Переменный размер шага, методы переменного порядка для уравнений с частными производными

ДеКария, Виктор (2019) Переменный размер шага, методы переменного порядка для уравнений с частными производными. Докторская диссертация, Питтсбургский университет. (Не опубликовано)

Это последняя версия данного товара.

Аннотация

Методы переменного размера шага и переменного порядка (VSVO) – это методы выбора для эффективного решения широкого диапазона ODE с минимальными трудозатратами и гарантированной точностью. Однако методы VSVO имеют ограниченное влияние в сложных приложениях из-за их вычислительной сложности и сложности их реализации в устаревшем коде. Целью данной диссертации является разработка, анализ и тестирование новых методов VSVO, которые имеют ту же вычислительную сложность, что и их неадаптивные аналоги на шаг.Адаптивность позволяет этим методам выполнять меньше действий, что делает их в целом менее сложными.

Здесь мы покажем, как использовать любой метод формулы обратного дифференцирования (BDF) в качестве основы для метода VSVO. Адаптивность к порядку достигается с помощью недорогого метода постобработки, известного как временная фильтрация. Временные фильтры не увеличивают асимптотическую сложность этих методов и позволяют вычислять каждый возможный порядок в семействе VSVO по той же цене, что и одно решение BDF. Этот подход дает новые нестандартные методы измерения временного шага, которых нет в литературе, и мы анализируем их стабильность и точность здесь.

Обратный Эйлер (BDF1) и BDF2 – чрезвычайно распространенные методы, и это исследование демонстрирует, как их можно преобразовать в последовательные адаптивные коды с помощью всего лишь нескольких дополнительных строк кода. Мы также разрабатываем решатель под названием Multiple Order One Solve Embedded 2,3,4 (MOOSE234). MOOSE234 – это метод VSVO, основанный на BDF3, который вычисляет приближения второго, третьего и четвертого порядка на каждом шаге. Все три приближения в MOOSE234 являются по крайней мере A (альфа) стабильными, а приближение второго порядка является стабильным.

Хотя эти методы обычно применимы к любой системе первого порядка по времени, мы сосредоточены на вопросах, относящихся к уравнениям Навье-Стокса. Наши методы были оптимизированы для решателей Навье-Стокса, и мы включаем линейно-неявные и неявно-явные (IMEX) версии.


Поделиться

Образец цитирования / Экспорт: Выберите формат … Цитата – TextCitation – HTMLEndnoteBibTexDublin CoreOpenURLMARC (ISO 2709) METSMODSEP3 XMLReference ManagerRefer
Социальные сети:

Детали

Тип изделия: Университет Питтсбурга ETD
Статус: Неопубликованные
Создатели / Авторы:
Комитет ETD:
Дата: 25 сентября 2019
Тип даты: Публикация
Дата защиты: 25 июля 2019
Дата утверждения: 25 сентября 2019
Дата отправки: 19 июля 2019
Ограничение доступа: Без ограничений; Немедленно выпустите ETD для доступа по всему миру.
Количество страниц: 110
Учреждение: Питтсбургский университет
Школы и программы: Школа искусств и наук Дитриха> Математика
Степень: PhD – доктор философских наук
Тип диссертации: Докторская диссертация
Реферировано: Есть
Неконтролируемые ключевые слова: Адаптивный временной шаг, обратный Эйлер, обратная формула дифференцирования, BDF, временной фильтр, переменный размер шага, переменный порядок
Дата депонирования: 25 сен 2019 14:56
Последнее изменение: 25 сен 2019 14:56
URI: http: // d-scholarship.pitt.edu/id/eprint/37251

Доступные версии этого товара


Метрики

просмотров по месяцам за последние 3 года

Plum Analytics


Действия (требуется логин)

Просмотреть товар

12. Численное решение Рунге-Кутта (RK4) для дифференциальных уравнений

В последнем разделе метод Эйлера дал нам один из возможных подходов к численному решению дифференциальных уравнений.

Проблема с методом Эйлера заключается в том, что для получения достаточно точного результата необходимо использовать небольшой размер интервала. То есть это не очень эффективно.

Метод Рунге-Кутта дает лучший результат за меньшее количество шагов.

Метод Рунге-Кутта Формула порядка 4

Заявки РК4

Механика
Биология
  • Хищник-жертва модели
  • Обвал рыболовства
  • Доставка лекарств
  • Предсказание эпидемии
Физика
  • Модели изменения климата
  • Озоновая защита
Авиация
  • Бортовые компьютеры
  • Аэродинамика

`y (x + h)` = y (x) + `1/6 (F_1 + 2F_2 + 2F_3 + F_4)`

где

`F_1 = hf (x, y)`

`F_2 = hf (x + h / 2, y + F_1 / 2)`

`F_3 = hf (x + h / 2, y + F_2 / 2)`

`F_4 = hf (x + h, y + F_3)`

Откуда взялась эта формула?

Вот краткая предыстория формулы.

Вывод формулы

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

Однако одним из недостатков метода Тейлора является то, что вам нужно дифференцировать функцию один раз для каждого нового члена, который вы хотите вычислить. Это может быть проблематично для сложных функций и плохо работает при компьютерном моделировании.

Карл Рунге (произносится как «рунга») и Вильгельм Кутта (произносится как «коота») стремились предоставить метод приближения функции без необходимости дифференцировать исходное уравнение.

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

Метод заказа 2 Рунге-Кутты

Начнем с двух вычислений функции вида:

`F_1 = hf (x, y)`

`F_2 = hf (x + alpha h, y + beta F_1)`

«Альфа» и «бета» – неизвестные величины. Идея заключалась в том, чтобы взять линейную комбинацию членов `F_1` и` F_2`, чтобы получить приближение для значения `y` при` x = x_0 + h` и найти подходящие значения `alpha` и` beta`. .

Путем сравнения значений, полученных с использованием метода рядов Тейлора и приведенных выше терминов (я избавлю вас от подробностей здесь), они получили следующее: Метод Рунге-Кутты Порядка 2 :

`y (x + h) = y (x) +1/2 (F_1 + F_2)`

где

`F_1 = hf (x, y)`

`F_2 = hf (x + h, y + F_1)`

Метод заказа Рунге-Кутты 3

Как обычно в этой работе, чем больше терминов мы возьмем, тем лучше будет решение.На практике решение порядка 2 используется редко, поскольку оно не очень точное.

Лучший результат дает метод Порядка 3:

`y (x + h) =` y (x) +1/9 (2F_1 + 3F_2 + 4F_3) `

где

`F_1 = hf (x, y)`

`F_2 = hf (x + h / 2, y + F_1 / 2)`

`F_3 = hf (x + (3h) / 4, y + (3F_2) / 4)`

Это было получено аналогично предыдущей формуле путем сравнения результатов ряда Тейлора.

Наиболее часто используемая формула Рунге-Кутта – это формула четвертого порядка (RK4), поскольку она дает лучший компромисс между вычислительными требованиями и точностью.(0,1 + 0,96545654732)) “ = -0,03154393258`

Шаг 2

Затем мы берем эти 4 значения и подставляем их в формулу Рунге-Кутты RK4:

`y (x + h)` = y (x) +1/6 (F_1 + 2F_2 + 2F_3 + F_4) `

`= 1 + 1/6 (-0.03678794411“ – \ 2xx0.03454223937“ – \ 2 xx0.03454345267“ {: – \ 0.03154393258) `

`= 0,9655827899`

Используя это новое значение `y`, мы начинаем заново, находя новые` F_1`, `F_2`,` F_3` и `F_4`, и подставляем их в формулу Рунге-Кутты.

Продолжаем этот процесс и составляем следующую таблицу значений Рунге-Кутты. (Я использовал электронную таблицу для получения таблицы. Использование калькулятора очень утомительно и подвержено ошибкам.)

24093983903567
x г `F_1 = h dy / dx`` x + h / 2` `y + F_1 / 2` `F_2` `y + F_2 / 2` `F_3` `x + h` `y + F_3` `F_4`
0 1 −0.0367879441 0,05 0,9816060279 -0,0345422394 0,9827288803 -0,0345434527 0,1 0,9654565473 -0,0315439326
0,1 0,9655827899 -0,0315443 0,15 0,9498106398 -0,0278769283 0,9516443257 −0.0278867954 0,2 0,9376959945 -0,023647342
0,2 0,937796275 -0,023648185 0,25 0,9259721824 -0,0189267761 0,9283328869 -0,0189548088 0,3 0,9188414662 -0,0138576597
0,3 0.9189181059 -0,0138588628 0,35 0,86745 -0,0084782396 0,89861 -0,0085314167 0,4 0,

66892

-0,0029773028
0,4 0,

21929
-0,0029786344 0,45 0,28756 0,0026604329 0. 0,002580704 0,5 0,28969 0,0082022376
0,5 0, 0,0082010354 0,55 0, 0,013727301 0,91992 34895 0,0136258867 0,6 0,9266857257 0,018973147
0.6 0,9267065986 0,0189722976 0,65 0,9361927474 0,02 40794197 0,9387463085 0,0239658709 0,7 0,9506724696 0,0287752146
0,7 0,9506796142 0,0287748718 0,75 0,9650670501 0,03 324 48616 0.967302045 0,0331305132 0,8 0,9838101274 0,0372312889
0,8 0,9838057659 0,0372315245 0,85 1,0024215282 0,0409408747 1,0042762033 0,0408359751 0,9 1.024641741 0,0441484563
0.9 1.024628046 0,0441492608 0,95 1.0467026764 0,0470593807 1.0481577363 0,0469712279 1 1.0715992739 0,04947
1 1.0715783953

Я не включил никаких значений после `y` в нижнюю строку, поскольку мы не будем их использовать.

Вот график найденных нами решений от «x = 0» до «x = 1».

Для интереса я расширил результат до `x = 10`, и вот получившийся график:

Упражнение

Решите следующее, используя RK4 (Метод Рунге-Кутты порядка 4) для `0 le x le 2`. Используйте размер шага h = 0,2:

`dy / dx = (x + y) sin xy`

`y (0) = 5`

Ответ

Шаг 1

У нас есть x_0 = 0 и y_0 = 5.

`F_1 = hf (x, y)` `= 0,2 ((0 + 5) sin (0) (5))` = 0`

Для `F_2` нам нужно знать:

`x + h / 2 = 0 + 0,2 / 2 = 0,1` и

`y + F_1 / 2 = 5 + 0/2 = 5`

Подставляем их в выражение `F_2`:

`F_2 = hf (x + h / 2, y + F_1 / 2)` = 0,2 ((0,1 + 5) sin (0,1) (5)) `= 0,484937`

Для `F_3` нам нужно знать:

`y + F_2 / 2` = 5 + 0,484937 / 2` = 5,24450702469`

Так

`F_3 = hf (x + h / 2, y + F_2 / 2)`

`= 0.2 ((0,1 + 5,24450702469) `{: xxsin (0,1) (5,24450702469))`

`= 0,535232`

Для `F_4` нам нужно знать:

`y + F_3“ = 5 + 0.535232` `= 5.535232`

Так

`F_4 = hf (x + h, y + F_3)`

`= 0,2 ((0,2 + 5,535232)` {: xxsin (0,2) (5,535232)) `

`= 1.02589

1`

Шаг 2

Затем мы берем эти 4 значения и подставляем их в формулу Рунге-Кутты RK4:

`y (x + h) = y (x) +1/6 (F_1 + 2F_2 + 2F_3 + F_4)`

`= 5 + 1/6 (0+` `2xx0.484937+ “ 2xx0.535232 ” {: + 1.02589

1) `

`= 5.5124008953`

Как и раньше, нам нужно взять это значение `y_1` и использовать новое значение` x_1 = 0,2`, чтобы найти следующее значение, `y_2`, и так далее до` x = 2`.

Ниже приводится таблица результирующих значений.

Вот таблица значений, которую мы получаем.

Посмотреть таблицу

x г `F_1 = h dy / dx`` x + h / 2` `y + F_1 / 2` `F_2` `y + F_2 / 2` `F_3` `x + h` `y + F_3` `F_4`
0 5 0 0.1 5 0,48494 5.2445070247 0,53523
0,2 5.53523
1.02589

0,2 5.5124008953 1.0194689011 0,3 6.0221353458 1,229424454 6,1271131223 1,2397613573 0,4 6.7521622525 0,6102193187
0,4 6.6070775356 0,670350862 0,5 6.9422529666 -0,4816655538 6.3662447588 -0,0570142602 0,6 6.5 500632755 -1,0142481364
0,6 6.3702013853 -0,8771352065 0.7 5,93 1633782 -1,1235642458 5.8084192624 -1,03691 0,8 5,3311976162 -1,1055304726
0,8 5,318

04

-1,0980514561 0,9 4,7698753724 -1,0356506358 4,80 10757825 -1,0539783626 1 4.2649227378 -0,9493143311
1 4,28 11304698 -0,9595184486 1,1 3,8013712455 -0,8453508288 3.8584550553 -0,8850171765 1,2 3,3961132933 -0,73891
1,2 3,42 1268202 -0,7592177155 1.3 3,0416593442 -0,6304549629 3.1060407205 -0,6882207337 1,4 2,7330474682 -0,5227646767
1,4 2,7680459044 -0,5581856458 1,5 2.4889530815 -0,4450766109 2.5455075989 -0,5066587622 1.6 2.2613871422 -0,354308929
1,6 2.2987183509 -0,3984550218 1,7 2.09949084 -0,3150804545 2,1411781236 -0,3672394563 1,8 1,9314788946 -0,2454079264
1,8 1.9639678892 -0.2886730166 1,9 1,8196313809 -0,230980612 1,8484775833 -0,2714612259 2 1.6925066634 -0,1779964411
2 1,71870

Еще раз, я не включил никаких значений после `y` в нижнюю строку, поскольку они не требуются.

Вот график найденных нами решений от «x = 0» до «x = 2».

Для интереса я расширил результат до `x = 6`, и вот получившийся график:

Заметим, что график не очень плавный. Если мы используем меньшее значение h, кривая в целом будет более гладкой.

Фактически, эта кривая является асимптотической по отношению к оси “x” (она плавно приближается к оси по мере увеличения “x”). На приведенном выше графике показано, что происходит, когда наши интервалы слишком грубые.(Этих движений после примерно `x> 4.5` не должно быть.)

Вот снова график решения, но на этот раз я использовал «h = 0,1». Это лучше, но все же есть “покачивание” около x = 5, которого не должно быть.

Если бы я взял еще меньшие значения `h`, я бы получил более плавную кривую. Однако это только до определенного момента, потому что со временем ошибки округления становятся значительными. Кроме того, время вычислений увеличивается без особой пользы.

Улучшения Рунге-Кутта

Как вы можете видеть из приведенного выше примера, на кривой есть точки, где значения `y` меняются относительно быстро (между` 0` и `1`), и другие места, где кривая более близка к линейной (от` 1` до `2`).Кроме того, у нас странное поведение около x = 5.

На практике, при написании компьютерной программы для выполнения Рунге-Кутты, мы допускаем переменные значения «h» – довольно маленькие, когда кривая изменяется быстро, и большие значения «h», когда кривая относительно гладкая.

Системы компьютерной алгебры Mathematics (такие как Mathcad, Mathematica и Maple) включают в себя процедуры, которые вычисляют RK4 за нас, поэтому нам просто нужно указать функцию и интересующий интервал.

Осторожно

Всегда будьте осторожны со своими ответами! Численные решения – это только приближения, и при определенных условиях вы можете получить хаотические результаты.

Нужна помощь в решении другой задачи с исчислением? Попробуйте решить проблемы.

Заявление об ограничении ответственности: IntMath.com не гарантирует точность результатов. Решение проблем, предоставленное Mathway.

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


Образование: Компьютерные науки и математика, Стэнфордский университет

Предметы: Алгебра, геометрия, триггер, исчисление, Adv.Исчисление, подготовка к SAT / ACT

Биография: Привет! Я Лигия, мне очень нравится наставлять других, и я работала репетитором математики в Стэнфорде. Я также интересуюсь физическими науками, и я провожу исследования с использованием машинного обучения в биологических проблемах. Я люблю теннис, бразильскую музыку и все виды рукоделия.

Стоимость: 100 $ / час

.

Оставить комментарий