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

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

Уфимский государственный авиационный технический университет

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

Курсовая работа

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

Уфа 2007 г.

Введение

Целью курсовой работы является: написать программу для нахождения приближенного решения обыкновенного дифференциального уравнения y'=f(x,y), y(a)=y0 методом Рунге-Кутта пятого порядка на отрезке [a,b] с заданным постоянным шагом h.

Для достижения данной цели необходимо выполнить следующие задачи:

1. Рассмотреть суть метода Рунге-Кутта.

2. Назначение и область применения.

3. Протестировать программу.

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

Для решения задачи будет использоваться язык программирования Turbo Pascal 7.0, так как этот язык позволяет работать с математическими формулами, проводить различного рода математические операции и действия. Turbo Pascal фирмы Borland является расширением стандарта языка и содержит интегрированную среду, данного ускоряющую и облегчающую процесс разработки программ. В языке программирования Turbo Pascal 7.0 используется типизированный адресный оператор, открытые массивы и строки, что предоставляет пользователю дополнительные возможности при решении математических задач. В математических задачах часто требуется реализовать численные методы, экспериментально исследовать условие и скорость сходимости методов. В условии задачи, как правило, дается основная идея каждого метода (Эйлера, Рунге-Кутта и т.д.).

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

Глава 1. Суть метода Рунге-Кутта

Метод Рунге-Кутта включает в себя несколько других таких как: метод Эйлера и метод Эйлера - Коши. Методы Рунге-Кутта обладают следующими свойствами:

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

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

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

Рассмотрим сначала геометрическое построение и выведем некоторые формулы на основе геометрических аналогий После этого мы подтвердим полученные результаты аналитически (Аналитический метод, применяется которого дает решение дифференциального уравнения в виде аналитического выражения; Графический метод, дающий приближенное решение в виде графика; Численный метод, когда искомая функция получается в виде таблицы.) Предположим нам известна точка xmym на искомой кривой Тогда мы можем провести прямую линию с тангенсом угла наклона уm=f(xmym) которая пройдет через точку xmym Это построение показано на рис1 где кривая представляет собой точное но конечно неизвестное решение уравнения а прямая линия L1 построена так как это только что описано



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

Уравнение прямой L1 выглядит так: y=ym+ym(x-xm) так как y=f(xmym) и кроме того xm+1=xm+h тогда уравнение примет вид

ym+1=ym+h*f(xmym) 1.1

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

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

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

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

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

Ф(xmymh)=[f(xmym)+f(xm+hym+ymh)] 1.2

Где, ym=f(xmym)

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

y=ym+(x-xm)Ф(xmymh)

ym+1=ym+hФ(xmymh) 1.3

Соотношения 12 13 14 описывают исправленный метод Эйлера



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

f(xy)=f(xmym)+(x-xm)f/x+(y-ym)f/x+ 1.5

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

Подставляя в формулу 15 x=xm+h и y=ym+hym и используя выражение 13 для ym получаем

f(xm+hym+hym)=f+hfx+hffy+O(h2)

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

Ф(xmymh)=f+h/2(fx+ffy)+O(h2)

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

ym+1=ym+hf+h2/2(fx+ffy)+O(h3)

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

Ф(xmymh)=f+(xm+h/2ym+h/2*ym) 1.6

где ym=f(xmym) 1.7

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

y=ym+(x-xm)Ф(xmymh)

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

ym+1=ym+hФ(xmymh)

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

ym+1=ym+hФ(xmymh) 1.9

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

Ф(xmymh)=a1f(xmym)+a2f(xm+b1hym+b2hym) 1.10

Где ym=f(xmym) 111

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

a1=a2=1/2; b1=b2=1

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

a1=0 a2=1

b1=b2=1/2

Формулы 19 110 111 описывают некоторый метод типа Рунге-Кутты Посмотрим какого порядка метод можно рассчитывать получить в лучшем случае и каковы допустимые значения параметров a1 a2 b1 и b2

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

В разложении f(xy) в ряд 15 в окрестности точки xmym положим x=xm+b1h y=ym+b2hf

Тогда f(xm+b1hym+b2hf)=f+b1hfx+b2hffy+O(h2) где функция и производные в правой части равенства вычислены в точке xmym

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

ym+1=ym+h[a1f+a2f+h(a2b1fx+a2b2ffy)]+O(h3)

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

ym+1=ym+h[a1f+a2f+h(a2b1fx+a2b2ffy)]+O(h3)

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

Сравнивая члены содержащие h2fx получаем a2b1=1/2

Сравнивая члены содержащие h2ffy получаем a2b2=1/2

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

Положим например a2=0 тогда a1=1- b1=b2=1/2 и соотношения 19 110 111 сведутся к

ym+1=ym+h[(1-)f(xmym)+f(xm+h/2ym+h/2f(xmym))]+O(h3) 1.12

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

et=kh3 1.13

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

метод рунге кутт программа уравнение

ym+1=ym+h/6(R1+2R2+2R3+R4) 1.14

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

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

1.1 Назначение и область применения

Данная задача может применяться в математических дисциплинах, таких как, численные методы, математические методы и т.д. Обширное семейство методов приближенного решения дифференциальных уравнений основано на ведении сетки и замене производных. Данный шаг позволяет избавиться от такого неудобного для компьютера объекта, как производная, заменив исходную задачу задачей алгебраической. В общем случае метод Рунге-Кутта является n-шаговым методом. На интервале [t0,t1] имеется последовательность n точек, и для этих точек вычисляются наклоны касательных к графику функции. Получив несколько значений, используем их взвешенное среднее для построения отрезка линии, идущего из начальной точки. Делается это для того, чтобы вычислить оценку решения в следующей точке. Затем вычисляется наклон в этой точке и т.д. Это - геометрическое истолкование методов Рунге-Кутта. Числовые коэффициенты выбираются из соображений точности; методы Рунге-Кутта различаются способом выбора этих коэффициентов. Методы Рунге - Кутта при m>5 не используются.

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

Найти приближенное решение обыкновенного дифференциального уравнения y'=f(x,y), y(a)=y0 методом Рунге-Кутта пятого порядка на отрезке [a,b] с заданным постоянным шагом h. Значение функции y(x) в узловых точках вычисляются по формуле:

Y [i+1]=y[1]+h/6(K[i]+2K[2]+2K[3]+K[4]), i=0,1,2,…

Где K[1]=f(x[i],y[i]);

K[2]=f(x[i]+h/2, y[i]+h/2K[1]);

K[3]=f(x[i]+h/2, y[i]+h/2K[2]);

K[4]=f(x[i]+h, y+hK[3]).

2.1 Выбор состава технических и программных средств

Для решения задачи «методом Рунге-Кутта» используется компьютер с минимальными требованиями:

1. Процессор - Intel 166.

2. ОЗУ 16Мб, видеокарта 8Мб.

3. Операционная система MS - DOS.

4. Винчестер (жесткий диск)-1Гб.

2.2 Вызов и загрузка программы

Данная программа находится на диске А(floppy disk).

Для запуска необходимо вставить floppy disk в дисковод.

На рабочем столе открываем «Мой компьютер», диск А, файл «Rkutt.exe»

2.3 Тестирование программы

Программа работает на языке Turbo Pascal 7.0, по методу Рунге-Кутта.

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

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

Данная программа приближённо вычисляет обыкновенные дифференцированные уравнения по методу Рунге-Кутта. Вводятся исходные данные (a,b,h1,k) и в результате появляются результативные данные (приложение №1).

Заключение

В ходе выполнения курсовой работы выполнена цель: найти приближенное решение обыкновенного дифференциального уравнения y'=f(x,y), y(a)=y0 методом Рунге-Кутта пятого порядка на отрезке [a,b] с заданным постоянным шагом h. Значение функции y(x) в узловых точках вычисляются по формуле:

Y [i+1]=y[1]+h/6(K[i]+2K[2]+2K[3]+K[4]), i=0,1,2,…

Где K[1]=f(x[i],y[i]);

K[2]=f(x[i]+h/2, y[i]+h/2K[1]);

K[3]=f(x[i]+h/2, y[i]+h/2K[2]);

K[4]=f(x[i]+h, y+hK[3]).

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

Список литературы

1. Бахвалов Н.С. Численные методы учебник для ВЗУов М., Наука, 1978

2. Воробьева Г.Н. Данилова А.Н. Практикум по численным методам: М., учебник 1987

3. Демидович Б.П. и Марон И.А. Основы вычислительной математики учебник для ВЗУов издание четвертое М., Наука, 1970

4. Каханер Дэвид. Численные методы и программное обеспечение учебник для ВЗУов М., Мир, 1998

5. Немнюгин С.А. Turbo Pascal.Практикум,учебник для ВУЗов, второе издание, С.-П., издание Питер. 2003.

6. Немнюгин С.А. Turbo Pascal. Программирование на языке высокого уровня. Учебник для ВУЗов, второе издание. С.-П., издание Питер, 2003.

7. Новикова Ф.А Дискретная математика для программистов. С.-П., издание Питер 2001.

8. Рюттен Т.Г. Франкен П.Р. Turbo Pascal 7.0 - К.: Торгово-издательское бюро BHV, 1996 - 448 С.: ил.

9. Турчак Л.И. Основы численных методов учебник для ВЗУов М., Знание, 1987

10. Фараонов В.В. Турбо Паскаль (в 3-х книгах). Книга 2. Библиотека Turbo Vision.-М.: Учебно-инжерненый центр МВТУ-ФЕСТО ДИДАКТИК, 1993.

Приложение

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

PROGRAM RKutt;

uses crt;

const

nmax=8;

type

vec=array[1..nmax] of Real;

var

h,a,b:Integer; i:Integer;

procedure der(x:Real; y:vec; var f:vec);

var

ki:array[1..5] of real;

begin

for i:=1 to 5 do

f[i]:=(h/6)*ki[i]+2*ki[2]+2*ki[3]+ki[4];

y[i+1]:=y[1]+h*y[i];

y[i]:=x*x+5;

end;

procedure RK4(n:Integer;x,h:Real; var y:vec);

var

i,j:Integer;

h1,h2,q:Real;

y0,y1,f:Vec;

begin

h1:=0;h2:=h1/2;

for i:=1 to n do begin

y0[i]:=y[i];

y1[i]:=y[i];

end;

for j:=1 to 4 do begin

der(x+h1,y,f);

if j=3 then h1:=h else h1:=h2;

for i:=1 to n do begin

q:=h1*f[i];

y[i]:=y0[i]+q;

y1[i]:=y1[i]+q/3;

end;

end;

for i:=1 to n do

y[i]:=y1[i];

end;

var

k,n:Integer;h1,y1:Real;

k1,k2,k3,k4:Integer;

x,ki:Array[1..5] of real;

y:Array[1..5] of real;

begin

writeln('Enter a=');

Readln(a);

writeln('Enter b=');

Readln(b);

writeln('Enter h1=');

Readln(h1);

writeln('Enter k=');

Readln(k);

n:=3;

h1:=a;

ki[1]:=k;

for i:=1 to n do begin

y[i+1]:=y1+(h1/6)*(Ki[i]+2*ki[2]+2*ki[3]+ki[4]);

writeln(y[i+1]);

end;

while h1<=b do begin

h1:=h1+h;

end;

Readln;

end.

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