Реферат: Решение дифференциальных уравнений. Обзор

--PAGE_BREAK--2 Решение дифференциальных уравнений с помощью Mathcad


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

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

набор точек в которых нужно найти решение;

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

Один из наиболее эффективных алгоритмов интегрирования ОДУ основан на численном методе Рунге-Кутты четвертого порядка. Функция, реализующая этот метод, имеет вид rkfixed (y,x1,x2, npoints,D)

Здесь:

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

x1, x2– граничные точки интервала, на котором ищется решение дифференциального уравнения. Начальные условия, заданные в векторе y,- это значение решения в точке x1;

npoints— число точек (не считая начальной точки), в которых ищется приближенное решение. При помощи этого аргумента определяется число строк (1+npoints) в матрице, возвращаемой функцией rkfixed;

D(x,y) – функция, возвращающая значение в виде вектора n элементов, содержащих первые производные неизвестных функций.
2.1 Метод Эйлера
Если задачу об отыскании всех решений дифференциального уравнения удается свести к конечному числу алгебраических операций, операций интегрирования и дифференцирования известных функций, то говорят, что уравнениеинтегрируется в квадратурах. В приложениях крайне редко встречаются уравнения, интегрируемые в квадратурах. Поэтому для исследования дифференциальных уравнений широко используются приближенные, численные методы их решения.

Численное решениена отрезке [a, b] задачи Коши
y' =f(x, y), y(a) =y0
состоит в построении таблицы приближенных значений
y,y1, ...,yi, ...yN
решенияy(x)в узлах сетки
a=x<x1<… <xi< ...<xN=b, y(xi)@yi.
Еслиxi=a+ i h, h=(b-a)/ N,то сетка <img width=«72» height=«27» src=«ref-1_1489603695-181.coolpic» alt=«Описание: t_1.gif (141 bytes)» v:shapes=«Рисунок_x0020_615»> называетсяравномерной.

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

Простейший одношаговый метод численного решения задачи Коши -метод Эйлера. В методе Эйлера величиныyiвычисляются по формуле



yi+1 =yi+hf(xi,yi), i= 0, 1
Найдем методом Эйлера на [0, 1] с шагом h=0.2 приближенное решение задачи Коши

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

Определим правую часть уравнения

Расчетные формулы метода Эйлера для решения этой задачи имеют вид
x0=0, y0= 1, xi+1 = xi + 0.2, yi+1 = yi + 0.2(sinxi — cosyi), i =0, 1, ..., 4.
Изобразим приближенное решение графически.
y' = sin x – cos y, y(0)=1.
Определим диапазон изменения номера точки i=0,1, ..., 4

Знак присваивания можно ввести щелчком по соответствующей позиции в панели Evaluation.
<img width=«146» height=«20» src=«ref-1_1489603876-1861.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex1/ex10135.JPG» v:shapes=«Рисунок_x0020_569»>
Определим начальное условие — решение в начальной точке

Для того чтобы ввести символ диапазона изменения индекса <..>, щелкните по соответствующей позиции в панели Matrix или введите с клавиатуры символ <;> («точка с запятой»)
<img width=«49» height=«20» src=«ref-1_1489605737-988.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex1/ex10142.JPG» v:shapes=«Рисунок_x0020_568»>
Определим шаг формулы Эйлера — шаг интегрирования

Для того чтобы ввести нижний индекс переменной, щелкните по соответствующей позиции в панели Matrix или в панели Calculator
<img width=«39» height=«21» src=«ref-1_1489606725-919.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex1/ex10217.JPG» v:shapes=«Рисунок_x0020_566»>

<img width=«39» height=«22» src=«ref-1_1489607644-931.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex1/ex10146.JPG» v:shapes=«Рисунок_x0020_567»>

Определим по формулам Эйлера значения приближенного решения в узлах сетки
<img width=«43» height=«20» src=«ref-1_1489608575-932.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex1/ex10193.JPG» v:shapes=«Рисунок_x0020_565»>
Выведем в рабочий документ вычисленные значения решения
<img width=«122» height=«23» src=«ref-1_1489609507-1640.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex1/ex10218.JPG» v:shapes=«Рисунок_x0020_563»>

<img width=«74» height=«23» src=«ref-1_1489611147-1120.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex1/ex10196.JPG» v:shapes=«Рисунок_x0020_564»>
Построим график найденного решения y(x)
<img width=«82» height=«133» src=«ref-1_1489612267-3210.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex1/ex10222.JPG» v:shapes=«Рисунок_x0020_561»>

<img width=«69» height=«133» src=«ref-1_1489615477-2474.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex1/ex10221.JPG» v:shapes=«Рисунок_x0020_562»>
Для того чтобы вывести значение переменной в рабочий документ, введите имя переменной, знак равенства и щелкните по рабочему документу вне выделяющей рамки

Для того чтобы построить график приближенного решения, щелкните в панели Graph по пиктограмме декартова графика, введите в помеченной позиции возле оси абсцисс обозначение компонент вектора, содержащего значения узлов сетки, а в позиции возле оси ординат — обозначение компонент вектора, содержащего значения приближенного решения в узлах сетки; затем щелкните по свободному месту в рабочем документе вне поля графиков.
<img width=«331» height=«163» src=«ref-1_1489617951-5504.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex1/ex10213.JPG» v:shapes=«Рисунок_x0020_559»>

<img width=«49» height=«20» src=«ref-1_1489623455-986.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex1/ex10224.JPG» v:shapes=«Рисунок_x0020_560»>

<img width=«336» height=«155» src=«ref-1_1489624441-7043.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex1/ex10226.JPG» v:shapes=«Рисунок_x0020_558»>
2.2 Метод Эйлера с шагом h/2.
Метод Эйлера допускает простуюгеометрическую интерпретацию. Пусть известна точка (xi,yi) интегральной кривой уравненияy'=f(x, y).

Касательная к интегральной кривой уравнения, проходящая через эту точку, определяется уравнением
y=yi+f(xi,yi)(x-xi).
Следовательно, вычисленная методом Эйлера точка (xi+1 ,yi+1 ),

Гдеxi+1=xi+h,yi+1=yi+h f(xi,yi), лежит на этой касательной.

Найдем методом Эйлера на [0, 1] с шагом h=0.2 и с шагом h=0.1 приближенное решение задачи Коши


y' = sin x – cosy,y(0)=1.
Изобразим приближенные решения графически.
Расчетные формулы метода Эйлера для решения этой задачи имеют вид
x0=0, y0= 1, xi+1 = xi + 0.2, yi+1 = yi + 0.2(sinxi — cosyi), i =0, 1, ..., 4

xi+1 = xi + 0.2, yi+1 = yi + 0.2(sinxi — cosyi), i =0, 1, ..., 9
Определим правую часть уравнения
<img width=«136» height=«19» src=«ref-1_1489631484-1833.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20038.JPG» v:shapes=«Рисунок_x0020_587»>
Знак присваивания можно ввести щелчком по соответствующей позиции в панели Evaluation.

Определим диапазон изменения номера точки i=0,1, ..., 4 для вычислений с шагом h=0.2
<img width=«45» height=«19» src=«ref-1_1489633317-997.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20041.JPG» v:shapes=«Рисунок_x0020_588»>
Для того чтобы ввести символ диапазона изменения индекса <..>, щелкните по соответствующей позиции в панели Matrix или введите с клавиатуры символ <;> («точка с запятой»)

При решении задачи с шагом h=0.2 назовем шаг h1, аргумент — x1, а решение — y1.

Определим начальное условие
<img width=«43» height=«21» src=«ref-1_1489634314-935.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20065.JPG» v:shapes=«Рисунок_x0020_589»>

<img width=«45» height=«21» src=«ref-1_1489635249-985.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20066.JPG» v:shapes=«Рисунок_x0020_590»>

Для того чтобы ввести нижний индекс переменной, щелкните по соответствующей позиции в панели Matrix или в панели Calculator

Определим шаг формулы Эйлера — шаг интегрирования
<img width=«49» height=«19» src=«ref-1_1489636234-1028.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20069.JPG» v:shapes=«Рисунок_x0020_591»>
Определим по формулам Эйлера значения приближенного решения в узлах сетки
<img width=«88» height=«21» src=«ref-1_1489637262-1395.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20071.JPG» v:shapes=«Рисунок_x0020_592»>

<img width=«147» height=«21» src=«ref-1_1489638657-2036.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20072.JPG» v:shapes=«Рисунок_x0020_593»>
Выведем в рабочий документ вычисленные значения решения

Для того чтобы вывести значение переменной в рабочий документ, введите имя переменной, знак равенства и щелкните по рабочему документу вне выделяющей рамки
<img width=«72» height=«128» src=«ref-1_1489640693-2727.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20150.JPG» v:shapes=«Рисунок_x0020_594»>

<img width=«85» height=«128» src=«ref-1_1489643420-3377.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20149.JPG» v:shapes=«Рисунок_x0020_595»>
Построим график найденного решения y1(x1)




<img width=«45» height=«19» src=«ref-1_1489646797-1017.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20151.JPG» v:shapes=«Рисунок_x0020_596»>

<img width=«317» height=«159» src=«ref-1_1489647814-5847.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20058.JPG» v:shapes=«Рисунок_x0020_597»>
Для того чтобы построить график приближенного решения, щелкните в панели Graph по пиктограмме декартова графика, введите в помеченной позиции возле оси абсцисс обозначение компонент вектора, содержащего значения узлов сетки, а в позиции возле оси ординат — обозначение компонент вектора, содержащего значения приближенного решения в узлах сетки; затем щелкните по свободному месту в рабочем документе вне поля графиков.

Определим диапазон изменения номера точки i=0,1, ..., 9 для вычислений с шагом h=0.1
<img width=«45» height=«19» src=«ref-1_1489653661-1031.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20095.JPG» v:shapes=«Рисунок_x0020_572»>
Для того чтобы ввести символ диапазона изменения индекса <..>, щелкните по соответствующей позиции в панели Matrix или введите с клавиатуры символ <;> («точка с запятой»)

При решении задачи с шагом h=0.1 назовем шаг h2, аргумент — x2, а решение — y2.

Определим начальное условие
<img width=«43» height=«21» src=«ref-1_1489654692-966.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20098.JPG» v:shapes=«Рисунок_x0020_573»>

<img width=«45» height=«21» src=«ref-1_1489655658-992.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20099.JPG» v:shapes=«Рисунок_x0020_574»>
Для того чтобы ввести нижний индекс переменной, щелкните по соответствующей позиции в панели Matrix или в панели Calculator

Определим шаг формулы Эйлера — шаг интегрирования
<img width=«49» height=«19» src=«ref-1_1489656650-1024.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20102.JPG» v:shapes=«Рисунок_x0020_575»>
Определим по формулам Эйлера значения приближенного решения в узлах сетки
<img width=«88» height=«21» src=«ref-1_1489657674-1435.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20104.JPG» v:shapes=«Рисунок_x0020_576»>

<img width=«147» height=«21» src=«ref-1_1489659109-2143.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20153.JPG» v:shapes=«Рисунок_x0020_577»>
Выведем в рабочий документ вычисленные значения решения. Для сравнения рядом выведены значения решения, вычисленные с большим шагом
<img width=«72» height=«128» src=«ref-1_1489640693-2727.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20160.JPG» v:shapes=«Рисунок_x0020_578»>

<img width=«85» height=«128» src=«ref-1_1489643420-3377.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20159.JPG» v:shapes=«Рисунок_x0020_579»>




<img width=«104» height=«242» src=«ref-1_1489667356-6810.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20132.JPG» v:shapes=«Рисунок_x0020_580»>
<img width=«121» height=«242» src=«ref-1_1489674166-8440.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20156.JPG» v:shapes=«Рисунок_x0020_581»>

Построим график решения y2(x2)
<img width=«51» height=«19» src=«ref-1_1489682606-1086.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20111.JPG» v:shapes=«Рисунок_x0020_582»>

<img width=«544» height=«175» src=«ref-1_1489683692-7148.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20112.JPG» v:shapes=«Рисунок_x0020_583»>

Построим на одном графике оба приближенные решения




<img width=«51» height=«19» src=«ref-1_1489682606-1086.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20161.JPG» v:shapes=«Рисунок_x0020_584»>

<img width=«531» height=«175» src=«ref-1_1489691926-8480.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20140.JPG» v:shapes=«Рисунок_x0020_585»>

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

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

Когда функции определены, щелкните по рабочему документу вне поля графиков.
<img width=«531» height=«175» src=«ref-1_1489700406-10690.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex2/ex20145.JPG» v:shapes=«Рисунок_x0020_586»>
2.3 Метод Рунге – Кутты
Методом Рунге-Кутты четвертого порядкаточности называют одношаговый метод, относящийся к широкому классу методов Рунге-Кутты. В этом методе величиныyi+1 вычисляются по следующим формулам:

yi+1 =yi+h(k1+ 2k2+ 2k3+k4)/6, i= 0, 1, ...

k1=f(xi,yi),

k2=f(xi+h/2,yi+hk1/2),

k3=f(xi+h/2,yi+hk2/2),

k4=f(xi+h,yi+hk3).
Найдем на [0, 1]приближенноерешение задачи Кошиy' = sinx– cosy,y(0)=1методом Рунге-Кутты 4-го порядка с шагом h=0.2 и методом Эйлера с тем же шагом.Изобразим оба приближенные решения графически

Для решения задачиметодом Рунге-Кутты воспользуемся функциейrkfixed

Определим начальное условие — решение в начальной точке
<img width=«39» height=«21» src=«ref-1_1489711096-920.coolpic» alt=«Описание: www.exponenta.ru/educat/class/courses/ode/theme2/mathcad/ex5/ex50091.JPG» v:shapes=«Рисунок_x0020_609»>    продолжение
--PAGE_BREAK--