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

Министерство образования и науки РФ

 

ПЕНЗЕНСКАЯ ГОСУДАРСТВЕННАЯ ТЕХНОЛОГИЧЕСКАЯ

АКАДЕМИЯ

 

Кафедра «Прикладная математика и исследование операций в экономике»

 

«УТВЕРЖДАЮ»

 

Зав. кафедрой ______________

 

«_____»______________20__ г.

 

 

 

 

 

 

 

КУРСОВАЯ РАБОТА

по дисциплине «Численные методы»

 

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

Вариант № 18

 

 

 

 

 

 

Автор работы:     Чичков Роман Геннадьевич

 

Специальность:                                       080116 («Математические методы в

                                                                                   экономике»)

Обозначение курсовой работы:            ПГТА 080116-08КР000.17 ПЗ

Группа:       08ММ

 

Руководитель:      Козлов А.Ю.

 

Работа защищена «__»_____20__г.  Оценка _______________

 

 

 

 

 

2011 г.

 

Содержание

 

 

 

Введение

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

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

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

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

Поставлена следующая задача:

1. Научиться решать дифференциальные уравнения численными и символьными методами с помощью программных средств Matlab.

 

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

1.1 Обыкновенные дифференциальные уравнения

Обыкновенные дифференциальные уравнения (ОДУ) — это дифференциальное уравнение вида F(t, x, x’ ,x’’, … , x(n))=0, где x(t) — неизвестная функция, зависящая от переменной t. Порядком уравнения называется максимальный порядок n входящей в него производной (или дифференциала).

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

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

Как следует из определения, обыкновенным дифференциальным уравнением первого порядка называется уравнение F(x, y,y’)=0, где x - независимая переменная, y(x) - неизвестная функция. В форме, разрешённой относительно производной, уравнение первого порядка записывается так: y’=f(x,y).

Общее решение (общий интеграл) уравнения при n = 1 имеет вид Ф(x,y,C)=0 или y=μ(x, C).

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

Задача Коши обычно возникает при анализе процессов, определяемых дифференциальным законом эволюции и начальным состоянием (математическим выражением которых и являются уравнение и начальное условие). Этим мотивируется терминология и выбор обозначений: начальные данные задаются при t = 0, а решение отыскивается при t > 0.

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

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

Система обыкновенных дифференциальных уравнений n –го порядка может быть записана в канонической форме:

в нормальной форме:

или в векторной форме:

Решением системы обыкновенных дифференциальных уравнений Y'=F(x,Y) называется вектор–функция Y(x)=Φ(x) , которая определена и непрерывно дифференцируема на промежутке (a; b) и удовлетворяет системе Y ' = F(x,Y) на этом промежутке.

Задачей Коши для системы обыкновенных дифференциальных уравнений называется следующая задача: найти решение Y(x) системы Y ' = F(x,Y) такое, что Y(x0) = Y0 . Здесь

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

Вектор–функция Y=Y(x, C)=Y(x, C1,C2, … , Cn) , зависящая от n произвольных постоянных C1,C2, … , Cn называется общим решением системы дифференциальных уравнений на [a; b], если:

при любых допустимых значениях постоянных C1,C2, … , Cn функция Y(x, C) является решением системы на [a;b];

какова бы ни была начальная точка (x0, Y0) из области определения правой части системы, существуют такие значения C*1,C*2, … , C*n постоянных C1, C2, … , Cn , что функция Y(x, C*1,C*2, … , C*n ) является решением задачи Коши Y(x0)=0.

Пусть Y(x)=Φ(x) — решение системы, определенное на [a, b]. Тогда множество точек {Φ(x)}, x Є[a, b] — кривая в пространстве Rn.

Эту кривую называют фазовой траекторией или просто траекторией системы, а пространство Rn , в котором расположены фазовые траектории, фазовым пространством системы.

2. Метод Рунге-Кутты

2.1 Описание метода

Существует множество методов решения дифференциальных уравнений и систем ОДУ, такие как метод Эйлера, метод Адамса, метод Нистрема, семейство методов Рунге-Кутты. На последнем остановимся попродробнее.

Методы Рунге - Кутты - важное семейство численных алгоритмов решения обыкновенных дифференциальных уравнений и их систем. Данные итеративные методы явного и неявного приближённого вычисления были разработаны около 1900 года немецкими математиками К. Рунге и М. В. Куттой.

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

Метод Рунге—Кутты 4 порядка столь широко распространён, что его часто называют просто методом Рунге—Кутта.

Рассмотрим задачу Коши

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

где h — величина шага сетки по x и вычисление нового значение проходит в четыре этапа:

Метод Рунге-Кутты имеет погрешность, пропорциональную h5 на каждом шаге численного интегрирования и поэтому является одним из наиболее употребительных методов численного решения задач Коши для ОДУ и систем ОДУ.

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

X’=f(t,X,Y,…),

Y’=f(t,X,Y,…)

и т.д.,

которые имеют решение:

Х=Х(t)

Y=Y(t)

и т.д.,

где t - независимая переменная (например, время); X, Y и т.д. - искомые функции (зависимые от t переменные). Функции f, g и т.д. - заданы. Также предполагаются заданными и начальные условия, т.е. значения искомых функций в начальный момент.

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

Метод Рунге-Кутта заключается в рекуррентном применении следующих формул:

,

, … ,

где

,

, … ,

, … ,

,

, … ,

,

, … 

2.2 Реализация метода в системе Matlab

Система MatLab имеет возможность решения систем обыкновенных дифференциальных уравнений. Для этого в ней имеются необходимые «решатели». Это функции ode23, ode45, odell3, odel5s, ode23s,ode23t и ode23tb.

Функции с суффиксом s предназначены для решения так называемых систем жестких дифференциальных уравнений. Для всех остальных систем дифференциальных уравнений наиболее употребительной является функция ode45, реализующая алгоритм Рунге- Кутта 4- 5-го порядка (разные порядки точности используются для контроля шага интегрирования).

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

yl'    = Fl( х,   yl,   у2,   ...,   уn );

у2'   = F2( х,   yl,   у2,   ...,   уn ) ;

уn'   = Fn( х,   yl,   у2,   ...,   уn);

то есть состоит из n уравнений, разрешенных относительно первых производных функций yl, у2, ..., уn (это все функции от х). Введем вектор-столбцы Y и F, состоящие из yl, у2, ..., уn и Fl, F2, ...,Fn соответственно. Тогда система дифференциальных уравнений примет следующий векторный вид:

Y' =F(х, Y);

Чтобы применить «решатель» ode45, нужно оформить в виде собственной функции правую часть системы уравнений F (х, Y).

 

3. Выполнение задания

Задание:

Решить систему дифференциальных уравнений методом Рунге - Кутты 4- 5-го порядка при заданных начальных и конечных условиях (постоянный h=0.1;0.01). Оценить точность результатов решения. Построить графики , , , а также график f(x,y,z) в декартовой системе координат. Создать видеофайл решения задачи.

№ п/п	СДУ	Начальные условия	Конечные условия
18.		0.5 0.0 0.5	20.0

 

 

Для данного примера неизвестная вектор-функция Y состоит из трех элементов:

Y = [yl, у2, y3]

Так как функции yl, у2 и y3 вычисляются во многих точках в процессе нахождения решения, то реально yl, у2 и y3 являются вектор-столбцами. Вектор F правых частей системы уравнений вычисляем с помощью собственной функции Div:

function F = div(t, y)

F(1) = 1./(y(1)^3.)*sin(y(2))-1./y(3),

F(2) = cos(sqrt(y(2)^2.+y(1)^3.))+2.*y(3),

F(3) = y(1)^2.*cos(2.*y(2)),

F = F';

end

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

Далее составим ещё одну функцию events для определения времени окончания выполнения программы

function [value,isterminal,direction] = events(t,Y)

global tf;

value = tf - t;

isterminal = 1;

direction = -1;

 

Теперь можно вызывать функцию ode45:

global tf;

tf = 2.9;

refine = 4;

options = odeset('Events',@events,'OutputFcn',@odeplot,'OutputSel',1,...

                 'Refine',refine);   

tspan = [0 tf];

  Y0 = [0.5;0;0.5];

[t,Y]  = ode45(@div, tspan,Y0, options);

t,

Y,

plot( t, Y(:,1));

hold on;

plot( t, Y(:,2),'r');

plot( t, Y(:,3),'g');

legend('x','y','z', 2);

grid on;

Решение отображается на графике (Рисунок 1.)

 

 

 

 

Рисунок 1 – Решение ОДУ

Данное решение получается при условии того, что h автоматически выбирается программой, исходя из поведения функции. Но шаг можно задать постоянный, и, независимо от поведения функции, он будет неизменным. Например, зададим шаг h=0,01. Текст программы в данном случае будет выглядеть так:

global tf;

 

 

tf = 0.8;

 

refine = 4;

options = odeset('Events',@events,'OutputFcn',@odeplot,'OutputSel',1,...

                 'Refine',refine);

 

            

tspan = [ 0: 0.01 : tf];

 

Y0 = [1;0;0];

[t,Y]  = ode45(@div, tspan,Y0, options);

t,

Y,

plot( t, Y(:,1));

hold on;

plot( t, Y(:,2),'r');

plot( t, Y(:,3),'g');

legend('x','y','z', 2);

grid on;

 

Получим иное решение (Рисунок 2.)

Рисунок 2 – Решение ОДУ с h=0.01

 

Различий по сравнению с первым графиком практически не наблюдается. Из-за того что h достаточно частый. Если же частоту шага уменьшить хотя бы до 0,1 то получим следующее: