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

Отчет по лабораторной работе

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

чиcленными методами

Как процесс математического

моделирования

Уравнения: 4

1. dy/dx = — (1 — y2)1/2 4

2. dy/dx = (y+(y2-4x2)1/2) / x 5

3. dy/dx = (4-y2)1/2 / y 6

4. dy/dx = 1-y4 7

5. dy/dx = 1/2×(y2-1) 8

6. x’=y-0.1x 9

y’=-x-0.1y 9

7. x’=y-x(x2+y2) 10

y’=-x-y(x2+y2) 10

8. dy/dx = xy(1+2y)/(1-y2) 11

Заключение. 12

Введение

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

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

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

Примечания:

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

Синий – метод Эйлера.

Зелёный – метод Штернера.

Фиолетовый – метод с регулировкой шага.

Белый – метод Рунге-Кутта

Красный – точное решение.

Уравнения: 1.  dy/dx = — (1 — y2)1/2

В данном примере при относительно больших шагах (H) происходит: нарушение достаточного условия единственности. Т. к. при Y стремящемся к ±1 происходит быстрая переориентация поля направлений, и касательная не успевает за производной функции f(x, y). Её можно избежать, уменьшив шаг (H).

Как это происходит?

На очередном шаге, очень близком к Yn → ±1 ломанная Эйлера не успевает отслеживать изгиб интегральной кривой, и вылетает за границу решения. На рисунке показано как это происходит.

Sn – имеющаяся ошибка, S(n+1) – ошибка на очередном шаге.

Постепенно уменьшая шаг интегрирования, я заметил, что метод Эйлера первый стал показывать безошибочное решение, затем метод Штермера, и только потом метод Рунге – Кутта. Метод с регулировкой шага в уравнении всегда дает ошибку, так как функция достаточно гладкая, этот метод с каждым новым расчетом увеличивает шаг и не может вовремя, в частности вблизи Y = ±1, «затормозить, из-за чего и вылетает за пределы существования функции.

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

Других видимых ошибок я в этом примере не обнаружил.

2. dy/dx = (y+(y2-4x2)1/2) / x

Рассмотрим нижнюю область. Близкую к Y=0.

Решая методом Эйлера и другими методами возни­кает ошибка : нарушение достаточных условий единст­венности.

т. е. быстрая пере­ориентация поля направле­ния. Причём всё это проис­ходит при достаточно боль­ших шагах.

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

Чем более Y → 0, тем меньшим должен быть шаг.

Теперь рассмотрим верхнюю область.

При H > 0 все методы без исключения дают правильный результат, так как направление поля очень ровное. Самое интересное начинается если H < 0. Здесь возможны варианты.

При малых шагах программа выдаёт результат: бесконечное зна­чение правой части: нарушение достаточного условия существования. Функция стремится к бесконечности, ошибка шага неограниченно возрастает.

При больших шагах возникает: нарушение достаточных условий единст­венности.

А если соответствующим образом подобрать шаг и начальную координату (см. рис.), то получим интересное явление. Функция проходит до Y=0, далее точно попадает в 0, тем самым она попадает в особую точку, и далее переходит на частное решение, по которой благополучно и уходит из поля нашего зрения. Функция получается словно «склеенной». «Шаг вправо, шаг влево» и вылезают вышеперечисленные ошибки.

Было замечено, что при x=0, y>0. программа вылетает.

3. dy/dx = (4-y2)1/2 / y

Сразу заметим, что прямые Y=±2 тоже удовлетво­ряют нашему уравнению, а Y = 0 не является ре­шением.

Решения полученные нашими методами совпадают с точным решением только когда Н > 0 ( — метод рег. шага).

Метод Эйлера выдаёт очень неточное решение. Метод Штермера и Рунге-Кутта близки к точному. С уменьшением шага ошибка нарушения достаточного условия единственности пропадает, и все ме­тоды начинают показывать достаточно точные решения ( — метод рег. шага).

Если H < 0, то происходит ошибка при переходе фазовых траекторий, связанная с нарушением достаточных условий существования.

Причём этим отличились все методы, и точного решения я так и не увидел ни при каких начальных условиях.

Так же было замечено, что чем ближе Y0 к Y=0, тем больший шаг доступен для правильного решения. Это можно объяснить тем, что с увеличением «длины» графика, касательной остается больше места для «стратегических маневров». Так же замечено, что чем выше порядок членов Н, тем большую разность шага можно допустить, это объясняется более точными вычислениями методов.

И всё это не относилось к методу регулировки шага.

4. dy/dx = 1-y4

При очень ма­лом шаге происхо­дит ошибка ко­нечно-разрядной арифметики. Ошибка происходит при ис­пользовании любых методов, кроме регулировки шага. При боль­ших шагах ошибок не происходит. Увеличив шаг до 0.00202, получаем достаточное совпа­дение. При уменьшении же шага ситуация ухудшается. Так же, чем ближе мы к Y = 1, тем большая ошибка обнаруживается, и чем далее от Y = 1, тем менее ее заметно (в окне графика функции). Это относиться к предельно маленьким шагам, недоступным к выставлению пользователем. При шаге 0.00202 такого эффекта уже не наблюдается.

При Y = 1, все методы правильные, причем при любом шаге ;).

Рассмотрим влияние ошибок конечно-разрядной арифметики при вы­полнении операции сложения:
предположим, что дифференциальное уравнение dy/dx = f (x, y) имеет частное решение Y=a¹0. Тогда это изоклина горизонтального наклона. Если она в то же время является асимптотой для решения некоторой задачи Коши с начальным условием Y(x0) = Y0¹a, то в процессе числен­ного интегрирования f(xn, yn)®0 при n®¥. Тогда h×f(xn, yn)®0, т. е., начиная с некоторого шага, это произведение станет очень мало, и при сложении с yn произойдет выравнивание порядков слагаемых; при этом значащий разряд этого произведения выйдет за пределы разрядной сетки и будет утерян. С этого момента к Yn будет прибавляться на каждом шаге нуль, и интегральная кривая выйдет на прямую, параллельную асимптоте, вместо того, чтобы неограниченно приближаться к ней. Заметим также, что с уменьшением h эта ошибка станет только больше.

5. dy/dx = 1/2×(y2-1)

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

X = -1,62566841

Y = 0,94285715

Шаг = 0,05

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

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

Других ошибок замечено не было.

6. x’=y-0.1x y’=-x-0.1y

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

При шаге h>0,198 вычислительная неустойчивость появляется и для метода Эйлера, что не скажешь про методы Штернера и Рунге-Кутта, которые продолжают считать правильно.

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

Рассмотрим процесс однажды возникшей ошибки. Пусть на шаге n ошибка равняется Jn. Тогда новое значение уn+1 будет подсчитано с ошибкой Jn+1.

yn+1 + Jn+1 =yn + jn + hf(xn, yn+jn) +o(h)

Считая jn, h малыми и отбрасывая нелинейные по члены в разложении Тейлора, в линейном приближении получим.

Вычислительная неустойчивость будет наблюдаться при таких h, при которых | 1+h ¶f(xn, yn)/¶y | принимает значения больше единицы. Для систем уравнений аналогичное условие заключается в том, что определитель матрицы

(Е +h Df(xn, yn)/Dy) (где Е – единичная матрица, а — матрица Якоби по переменным у) должен стать больше единицы.

Для данной системы Р находится так: Её определитель:

Условие | | >1 будет выполнено при h> 0.2/1.01 » 0.2

7. x’=y-x(x2+y2) y’=-x-y(x2+y2)

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

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

8. dy/dx = xy(1+2y)/(1-y2)

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

1) Y<-1, X>0, H>0

2) Y<-1, X<0, H<0

3) -1<Y<0 , любой X

если h>0 график вправо

если h<0 график влево

4) Y>1, X ≥ 0.60962558, h<0

5) Y>1, X ≤ 0.60962558, h<0

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

Выше прямой y=1, можно утверждать, что все численные методы допускают грубую ошибку при «y» близком к 1. т. к. тут нарушается достаточное условие существования решения.

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

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

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

Заключение.

В результате проведенных исследований можно подвести итоги проделанной работы:

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

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

— Научились ликвидировать эти ошибки различными методами.

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

— Узнали, что метод Рунге-Кутта наиболее точный, но медленный, а метод Эйлера наиболее быстрый, но менее точный. Хотя при современных технологиях разницы в выполнении расчетов я не заметил. Но, но…

…было интересно %)