Размещено на http:www.allbest.ru/

Федеральное агентство по образованию

Государственное образовательное учреждение

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

им. И.И. Ползунова»

БИЙСКИЙ ТЕХНОЛОГИЧЕСКИЙ ИНСТИТУТ (филиал)

Кафедра информатики и вычислительной математики

УДК 681.3.06

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

по дисциплине «Информатика»

Бийск 2010

Содержание

Введение

1 Теоретическая часть

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

1.2 Метод Рунге-Кутта

2 Практическая часть

2.1 Задание

2.2 Алгоритм решения

2.3 Листинг программы

2.4 Результат работы программы

Заключение

• Список использованных источников

Введение

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

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

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

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

Задание предполагает:

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

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

1 Теоретическая часть

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

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

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

( x, y, y¹, ... y⁽ⁿ⁾ )=0. 1.1

Из теории ОДУ известно, что уравнение (1.1) эквивалентно системе n уравнений первого порядка

_k(x, y_1, y₁' ,y₂ ,y₂ ', ... ,y_n ,y_n ')=0. 1.2

где k=1, ... , n.

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

Первый тип - это задачи Коши, или задачи с начальными условиями. Для таких задач кроме исходного уравнения (1.1) в некоторой точке x₀ должны быть заданы начальные условия, т.е. значения функции y(x) и ее производных

y(x₀)=y₀' , y'(x₀)=y₁₀, ... , y^(n-1)(x₀)=y_n-1,₀.

Для системы ОДУ типа (1.2) начальные условия задаются в виде

y₁(x₀)=y₁₀ , y₂(x₀)=y₂₀, ... , y_n(x₀)=y_n0. 1.3

Ко второму типу задач относятся так называемые граничные, или краевые задачи, в которых дополнительные условия задаются в виде функциональных соотношений между искомыми решениями. Количество условий должно совпадать с порядком n уравнения или системы. Если решение задачи определяется в интервале x є x₀ ,x_k, то такие условия могут быть заданы как на границах, так и внутри интервала. Минимальный порядок ОДУ, для которых может быть сформулирована граничная задача, равен двум.

Третий тип задач для ОДУ - это задачи на собственные значения. Такие задачи отличаются тем, что кроме искомых функций y(x) и их производных в уравнения входят дополнительно m неизвестных параметров ₁₂х_m которые называются собственными значениями Для единственности решения на интервале [x₀x_k] необходимо задать m+n граничных условий В качестве примера можно назвать задачи определения собственных частот коэффициентов диссипации структуры электромагнитных полей и механических напряжений в колебательных системах задачи нахождения фазовых коэффициентов коэффициентов затухания распределения напряженностей полей волновых процессов и тд

К численному решению ОДУ приходится обращаться когда не удается построить аналитическое решение задачи через известные функции Хотя для некоторых задач численные методы оказываются более эффективными даже при наличии аналитических решений

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

1.2 Метод Рунге-Кутта

Все методы Рунге-Кутта являются одношаговыми методами: для получения решения y_i+1 при x_i+1 надо знать решение в одной точке (y_i, x_i) . Это позволяет начать решение используя начальные условия. Указанная особенность допускает изменение шага интегрирования в любой точке в процессе счета, что дает возможность строить численные алгоритмы с автоматическим выбором шага.

Методы Рунге-Кутта обладают следующими свойствами:

1 Эти методы являются одноступенчатыми: чтобы найти у_m+1 нужна информация о предыдущей точке x_my_m

2 Они согласуются с рядом Тейлора вплоть до членов порядка hp где степень р различна для различных методов и называется порядковым номером или порядком метода

3 Они не требуют вычисления производных от f(xy) а требуют вычисления самой функции

Рассмотрим сначала геометрическое построение и выведем некоторые формулы на основе геометрических аналогий После этого мы подтвердим полученные результаты аналитически

Предположим нам известна точка x_my_m на искомой кривой Тогда мы можем провести прямую линию с тангенсом угла наклона у_m=f(x_my_m) которая пройдет через точку x_my_m Это построение показано на рис1 где кривая представляет собой точное но конечно неизвестное решение уравнения а прямая линия L₁ построена так как это только что описано



Рис. 1

Тогда следующей точкой решения можно считать ту где прямая L₁ пересечет ординату проведенную через точку x=x_m+1=x_m+h

Уравнение прямой L₁ выглядит так: y=y_m+y_m(x-x_m) так как y=f(x_my_m) и кроме того x_m+1=x_m+h тогда уравнение примет вид

y_m+1=y_m+h*f(x_my_m) 14

Ошибка при x=x_m+1 показана в виде отрезка е Очевидно найденное таким образом приближенное значение согласуется с разложением в ряд Тейлора вплоть до членов порядка h так что ошибка ограничения равна

е _t=Кh2

Хотя точка на рис. 1 была показана на кривой в действительности y_m является приближенным значением и не лежит точно на кривой

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

Рассмотрим исправленный метод Эйлера и модификационный метод Эйлера В исправленном методе Эйлера мы находим средний тангенс угла наклона касательной для двух точек: x_my_m и x_m+hy_m+hy_m Последняя точка есть та самая которая в методе Эйлера обозначалась x_m+1y_m+1 Геометрический процесс нахождения точки x_m+1y_m+1 можно проследить по рис2 С помощью метода Эйлера находится точка x_m+hy_m+hy_m лежащая на прямой L₁ В этой точке снова вычисляется тангенс дает прямую Наконец через точку x_my_m мы проводим прямую L параллельную Точка в которой прямая L пересечется с ординатой восстановленной из x=x_m+1=x_m+h и будет искомой точкой x_m+1y_m+1

Тангенс угла наклона прямой и прямой L равен

Ф(x_my_mh)=[f(x_my_m)+f(x_m+hy_m+y_mh)] 15

где y_m=f(x_my_m) 16

Уравнение линии L при этом записывается в виде

y=y_m+(x-x_m)Ф(x_my_mh)

так что

y_m+1=y_m+hФ(x_my_mh) 17

Соотношения 15 16 17 описывают исправленный метод Эйлера



Рис. 2

Чтобы выяснить насколько хорошо этот метод согласуется с разложением в ряд Тейлора вспомним что разложение в ряд функции f(xy) можно записать следующим образом:

f(xy)=f(x_my_m)+(x-x_m)f/x+(y-y_m)f/x+ 18

где частные производные вычисляются при x=x_m и y=y_m

Подставляя в формулу 18 x=x_m+h и y=y_m+hy_m и используя выражение 16 для y_m получаем

f(x_m+hy_m+hy_m)=f+hf_x+hf_y+O(h2)

где снова функция f и ее производные вычисляются в точке x_my_m Подставляя результат в 15 и производя необходимые преобразования получаем

Ф(x_my_mh)=f+h/2(f_x+f_y)+O(h2)

Подставим полученное выражение в 17 и сравним с рядом Тейлора

y_m+1=y_m+hf+h2/2(f_x+f_y)+O(h3)

Как видим исправленный метод Эйлера согласуется с разложением в ряд Тейлора вплоть до членов степени h2 являясь таким образом методом Рунге-Кутты второго порядка

Рассмотрим модификационный метод Эйлера Рассмотрим рис3 где первоначальное построение сделано так же как и на рис2 Но на этот раз мы берем точку лежащую на пересечении этой прямой и ординатой x=x+h/2 На рисунке эта точка образована через Р а ее ордината равна y=y_m+(h/2)y_m Вычислим тангенс угла наклона касательной в этой точке

Ф(x_my_mh)=f+(x_m+h/2y_m+h/2*y_m) 19

где y_m=f(x_my_m) 110

Прямая с таким наклоном проходящая через Р обозначена через * Вслед за тем мы проводим через точку xmym прямую параллельную * и обозначаем ее через L₀ Пересечение этой прямой с ординатой x=x_m+h и даст искомую точку x_m+1y_m+1 Уравнение прямой можно записать в виде

y=y_m+(x-x_m)Ф(x_my_mh)

где Ф задается формулой 16 Поэтому

y_m+1=y_m+hФ(x_my_mh) 111

Соотношения 19 110 111 описывают так называемый модификационный метод Эйлера и является еще одним методом Рунге-Кутта второго порядка Обобщим оба метода Заметим что оба метода описываются формулами вида

y_m+1=y_m+hФ(x_my_mh) 112

и в обоих случаях Ф имеет вид

Ф(x_my_mh)=a₁f(x_my_m)+a₂f(x_m+b₁hy_m+b₂hy_m) 113

где y_m=f(x_my_m) 114

В частности для исправленного метода Эйлера

a₁=a₂=1/2;

b₁=b₂=1

Рис. 3

В то время как для модификационного метода Эйлера

a₁=0 a₂=1

b₁=b₂=1/2

Формулы 112 113 114 описывают некоторый метод типа Рунге-Кутты Посмотрим какого порядка метод можно рассчитывать получить в лучшем случае и каковы допустимые значения параметров a₁ a₂ b₁ и b₂

Чтобы получить соответствие ряду Тейлора вплоть до членов степени h в общем случае достаточно одного параметра Чтобы получить согласование вплоть до членов степени h2 потребуется еще два параметра так как необходимо учитывать члены h2f_x и h2f_y Так как у нас имеется всего четыре параметра три из которых потребуются для создания согласования с рядом Тейлора вплоть до членов порядка h2 то самое лучшее на что здесь можно рассчитывать - это метод второго порядка

В разложении f(xy) в ряд 18 в окрестности точки x_my_m положим

x=x_m+b₁h

y=y_m+b₂hf

Тогда f(x_m+b₁hy_m+b₂hf)=f+b₁hf_x+b₂hf_y+O(h2) где функция и производные в правой части равенства вычислены в точке x_my_m

Тогда 112 можно переписать в виде

y_m+1=y_m+h[a₁f+a₂f+h(a₂b₁f_x+a₂b₂f_y)]+O(h³)

Сравнив эту формулу с разложением в ряд Тейлора можно переписать в виде

y_m+1=y_m+h[a₁f+a₂f+h(a₂b₁f_x+a₂b₂ff_y)]+O(h3)

Если потребовать совпадения членов hf то a₁+a₂=1

Сравнивая члены содержащие h2f_x получаем a₂b₁=1/2

Сравнивая члены содержащие h2f_y получаем a₂b₂=1/2

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

Положим например a₂=0 тогда a₁=1- b₁=b₂=1/2 и соотношения 112 113 114 сведутся к

y_m+1=y_m+h[(1-)f(x_my_m)+f(x_m+h/2y_m+h/2f(x_my_m))]+O(h3) 115

Это наиболее общая форма записи метода Рунге-Кутта второго порядка При =1/2 мы получаем исправленный метод Эйлера при =1 получаем модификационный метод Эйлера Для всех отличных от нуля ошибка ограничения равна

e_t=kh3 116

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

y_m+1=y_m+h/6(R₁+2R₂+2R₃+R₄) 117

где R₁=f(x_my_m) 118

R₂=f(x_m+h/2y_m+hR₁/2) 119

R₃=f(x_m+h/2y_m+hR₂/2) 120

R₄=f(x_m+h/2y_m+hR₃/2). 121

Ошибка ограничения для этого метода равна e_t=kh5

так что формулы 117-121 описывают метод четвертого порядка Заметим что при использовании этого метода функцию необходимо вычислять четыре раза

2 Практическая часть

2.1 Задание

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

Исходные данные:

Дифференциальное уравнения первого порядка:

,

где T=2; k=1; x(t)=1.

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

.

2.2 Алгоритм решения

Основная программа

Размещено на http:www.allbest.ru/

Вычисление функции

Размещено на http:www.allbest.ru/

Алгоритм Рунге-Кутта

Размещено на http:www.allbest.ru/

2.3 Листинг программы

Uses Crt;

const

k=1;

T=2;

type

TFunc = function(x,y: extended): extended;

Var

f1: Text;

n: word; {количество точек разбиения}

x0,x_end,y0: extended;

function func(x,y: extended): extended; far;

begin

func := (k*x-y)/T;

end;

function RungeKutt(f: TFunc; x0_,x_end,y0_: extended; n: word): extended;

{ где x0_ и y0_ - начальное условие

x_end - точка, в которой необходимо вычислить результат

n - количество шагов для вычисления результата }

var

i : word; { счетчик цикла }

x,h : extended; { текущая точка и длина шага }

res : extended; {переменная для накопления конечного результата функции }

k1,k2,k3,k4: extended; {вспомогательные переменные вычисления результата }

begin

h:= (x_end - x0_)/n; { Находим длину шага }

res:= y0_; { устанавливаем начальные значения}

x:=x0_;

writeln(F1,x:0:3,' ',res:0:3); {Вывод расчитанных значений в файл}

for i:=1 to n do

begin { вычисляем результат по методу Рунге-Кутта 4го порядка }

k1:=f(x,res);

k2:=f(x+h/2,res+h*k1/2);

k3:=f(x+h/2,res+h*k2/2);

k4:=f(x+h,res+h*k3);

res:=res+h*(k1+2*k2+2*k3+k4)/6;

x:=x+h; { переходим к следующей точке }

writeln(F1,x:0:3,' ',res:0:3); {Вывод расчитанных значений в файл}

end;

RungeKutt:=res; { присваиваем конечный результат функции }

end;

begin

ClrScr;

assign(F1,'C:\runge.txt'); {Настройка записи в текстовый файл}

rewrite(F1);

writeln('Численное решение дифференциальных уравнений:');

writeln('Алгоритм Рунге-Кутты 4-го порядка');

writeln(' Tdy(y)/dt + y(t)=kx T=2; k=1');

writeln;

{Считывание данных для расчета}

write('Введите x0='); readln(x0);

write('Введите x_end='); readln(x_end);

write('Введите n='); readln(n);

write('Введите y0='); readln(y0);

writeln;

writeln(' x y') {Вывод "шапки" таблицы на экран}

writeln(F1,' x y'); {Вывод "шапки" таблицы в фай}

writeln('n=',n, ' Результат: ',RungeKutt(func, x0, x_end, y0, n):5:5);

Close(F1);

readln;

end.

Заключение

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

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

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

Список использованных источников

Самарский А.А., Гулин А.В. Численные методы - М.: Наука, ГРФМЛ, 1989.- 432с.

Бахвалов Н.С. Численные методы - М.: Наука, ГРФМЛ,1987.-600 с.

Маликов В.Т. Вычислительные методы и применение - Киев: В.Ш.-1989.-213 с.

1. Мудров А.Е.Численные методы для ПЭВМ на языках Паскаль, Фортран и Бейсик - Томск, МП “Раско”, 1991 г.

Размещено на Allbest.ru

...