67723

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

Курсовая

Информатика, кибернетика и программирование

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

Русский

2016-08-04

257.5 KB

12 чел.

ФЕДЕРАЛЬНОЕ АГЕНСТВО СВЯЗИ

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ ОБРАЗОВАТЕЛЬНОЕ БЮДЖЕТНОЕ УЧРЕЖДЕНИЕ

ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ

"СИБИРСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ ТЕЛЕКОММУНИКАЦИЙ И

ИНФОРМАТИКИ"

УРАЛЬСКИЙ ТЕХНИЧЕСКИЙ ИНСТИТУТ СВЯЗИ И ИНФОРМАТИКИ(филиал)

УрТИСИ ФГОБУ ВПО "СибГУТИ"

КАФЕДРА "ВЫСШАЯ МАТЕМАТИКА И ФИЗИКА "

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

по информатике:

«Визуализация численных методов.

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

Вариант № 1

                                                                               Выполнил: Алексеев А.А.

                                                                                                               гр. ИТЕ-12Б

                                                                                    Проверил: Бикбулатова Н.Г.

Екатеринбург

2012 г.

Содержание:

Введение………………………………………………………………….3

1. Постановка задачи…………………………………………………….4

2. Описание методов решения…………………………………………..5

2. 1. Численные методы решения задачи Коши……………………….5

2. 2. Метод Эйлера ……………………..…………………………….5

2. 3. Метод Рунге-Кутта 4-го порядка………………………………….7

2. 4. Решение поставленной задачи методами Эйлера и

Рунге-Кутта 4-го порядка …………………………………….8

2. 4. 1. Метод Эйлера ………………….. ……….……………………8

2. 4. 2. Метод Рунге-Кутта 4-го порядка ……………………………9

3. Алгоритм решения задачи…………………………………………...10

3. 1. Алгоритмы подпрограмм.………………………………………....10

3. 1. 1. Подпрограмма метода Эйлера……………………..………..10

3. 1. 2. Подпрограмма метода Рунге-Кутта 4-го порядка ………..11

3. 1. 3. Подпрограмма общего решения ………..…………………12

3. 2. Алгоритм функции…………………………………………………12

3. 3. Алгоритм программы………………………………………………13

4. Форма программы…………………………………………………….16

5. Листинг программы…………………………………………………..17

6. Решение задачи в MathCad…………………………………………..19

Заключение………………………………………………………………21


                                                    
Введение

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

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

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

  1.  Написать программу для решения данного дифференциального уравнения двумя численными методами в программе Visual Basic.
  2.  Проверить решение с помощью приложения MathCad.
  3.  Сравнить полученные разными методами результаты с общим решением.

  1.  
    Постановка задачи

Решить методами Эйлера и Эйлера модифицированного задачу Коши для дифференциального уравнения 1-го порядка на отрезке [X0; Xk] с шагом h и начальным условием: Y(X0) = Y0.

Ответ должен быть получен в виде таблицы результатов:

X

Y(1)

Y(2)

YT

X0

Y0(1)

Y0(2)

Y(X0)

X1

Y1(1)

Y1(2)

Y(X1)

Xk

Yk(1)

Yk(2)

Y(Xk)

Где Y(1), Y(2) – решения, полученные различными численными методами, YT – точное решение дифференциального уравнения.

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

Данные таблицы визуализировать на форме в виде графиков.

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

  

Дифференциальное уравнение

X0

Xk

h

Y0

Общее решение

y’ = - x·y/(x+1)

1,2

2

0,1

1

y=c·(x+1) ·exp(-x)


2. Описание методов решения

2. 1. Численные методы решения задачи Коши

При использовании численных методов выполняется замена отрезка [х0, X] - области непрерывного изменения аргумента х множеством . состоящего из конечного числа точек х0 < х1 < ... < xn = Х - сеткой.

При этом xi называют узлами сетки.

Во многих методах используются равномерные сетки с шагом:

Задача Коши, определённая ранее на непрерывном отрезке [х0, X], заменяется её дискретным аналогом - системой уравнений, решая которую можно последовательно найти значения y1, y2,…,yn - приближённые значения функции в узлах сетки.

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

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

Методы прогноза и коррекции (многошаговые), в которых для отыскания следующей точки кривой у = f(x) требуется информация более чем об одной из  предыдущих точек.   Чтобы  получить достаточно точное  численное значение, часто прибегают к итерации. К числу таких методов относятся методы Милна, Адамса - Башфорта и Хемминга.

Явные методы, в которых функция Ф не зависит от yn+1.

Неявные методы, в которых функция Ф зависит от yn+1.

2.2 Метод Эйлера.

Иногда  этот  метод  называют   методом  Рунге-Кутта  первого   порядка точности.

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

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

Y’ = f(x, y)

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

y(x0) = y0

Выберем шаг h и введём обозначения:

xi = х0 + ih  и yi = y(xi),   где   i = 0, 1, 2, ...,

xi - узлы сетки,

yi - значение интегральной функции в узлах.

Иллюстрации к решению приведены на рисунке 2.

Проведем прямую АВ через точку (xi, yi) под углом α. При этом tg α = f(xi, yi)

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

Тогда yi+1 = yi + Δy

Из прямоугольного треугольника ABC  

Приравняем правые части tg α = f(xi, yi) и . Получим

Отсюда Δу = hf(xi, yi).

Подставим в это выражение формулу yi+1 = yi + Δy, а затем преобразуем его. В результате получаем формулу расчета очередной точки интегральной функции:

.

Рисунок 2. Метод Эйлера.

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

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

F(x, у) - заданная функция – должна

быть описана отдельно.

Входные параметры:

Х0, XK—начальное и конечное

значения независимой переменной;

Y0 – значение y0 из начального условия

y(x0) = y0;

N - количество отрезков разбиения;

Выходные параметры:

У - массив значений искомого решения

в узлах сетки.

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

Метод Эйлера - один из простейших методов численного решения обыкновенных дифференциальных уравнений. Но существенным его недостатком является большая погрешность вычислений. На рисунке 2 погрешность вычислений для io шага обозначена ε. С каждым шагом погрешность вычислений увеличивается.

2.3 Метод Рунге-Кутта 4 порядка

Пусть дано дифференциальное уравнение первого порядка   с начальным условием y(x0)=y0.   Выберем шаг h и введем обозначения:

xi = x0 + ih  и yi = y(xi),   где   i = 0, 1, 2, ... .

Аналогично описанному выше методу производится решение

дифференциального уравнения. Отличие состоит в делении шага на 4 части.

         Согласно методу Рунге-Кутта четвертого порядка, последовательные значения yi   искомой функции y определяются по формуле:

yi+1 = yi +∆yi                где i = 0, 1, 2 ...

             ∆y=(k1+2*k2+2*k3+k4)/6                           

                    

a числа k1, k2 ,k3, k4 на каждом шаге вычисляются по формулам:

k1  = h*f(xi, yi )

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

k3  = F(xi +h/2, yi +k2 /2)*h

k4  = F(xi +h, yi +k3 )*h

Это явный четырехэтапный метод 4 порядка точности.

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

F(x, у) - заданная функция - должна

быть описана отдельно.

Входные параметры:
Х0,
XК - начальное и конечное

значения независимой

переменной;

Y0 – значение y0 из начального условия

y(x0)=y0;

N - количество отрезков разбиения;

Выходные параметры:

Y - массив значений искомого решения

в узлах сетки.

2. 4. Решение поставленной задачи методами Эйлера и Рунге-Кутта 4 порядка

2. 4. 1. Метод Эйлера

1. Строим оси координат;

2. Отмечаем A(1,2; 1) – первую точку интегральной кривой;

3. Ищем угол наклона касательной к графику в точке A:

4. Строим касательную l0 в точке А под углом α0;

5. Находим х1 по формуле: xi = х0 + ih, где h – шаг интегрирования

x1 =1,2 + 1 · 0,1 = 1,3

6. Проводим прямую x = x1 = 0,1  до пересечения с прямой l0, отмечаем точку B(x1; y1);

7. Ищем y точки B:

Из прямоугольного треугольника ABC ,

Δy = y1 y0,

Δx = x1x0 = h,

f(x0; y0) = (y1y0)/h =>

y1 = y0 + h · (f(x0; y0)) = 1 + 0,1 · f(1,2;1) = 1-0,545454 = 0,945

Следовательно, точка B имеет координаты (1,3; 0,945).

2.4.2. Метод Рунге-Кутта 4 порядка

1. Строим оси координат;

2. Отмечаем А(1,2; 1) – первую точку интегральной кривой;

3. Ищем угол наклона касательной к графику в точке A:

4. Строим касательную l0 в точке А под углом α0;

5. Находим х1 по формуле: xi = х0 + ih

x1 = 1,2 + 1 · 0,1 = 1,3;

  1.  Находим по формулам:

k1=0,1·f(1,2;1)=0,1·(-0,545454)= - 0,0545

k2=0,1· f(1,2+0,1/2;1-0,0545/2)= -0,0540

k3=0,1· f(1,2+0,1/2;1-0,0540/2)= - 0,0540

k4=0,1· f(1,2+0,1;1-0,0540)= - 0,0535

y1=(- 0,0545+2·(-0,0540)+2·(-0,0540) - 0,0535)/6= - 0,054

 y2=1- 0,054=0,946

Следовательно, следующая точка графика решения имеет координаты (1,3; 0,946)


3. Алгоритм решения задачи

3. 1. Алгоритмы подпрограмм

3.1.1 Подпрограмма метода Эйлера

3.1.2 Подпрограмма метода Рунге-Кутта 4 порядка

3. 1. 3. Подпрограмма общего решения

3. 2. Алгоритм функции


3. 3.
Алгоритм программы

  1.   Форма программы

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

1111

Dim j() As Single

Dim x() As Single

Dim y() As Single

Dim o() As Single

Private n, i As Integer

Private xk, x0, kx, ky As Single

Private k, k1, k2, k3, k4 As Single

Private h, max, min, y0 As Single

Private Function f(a, b As Single) As Single

f = -a * b / (a + 1)

End Function

Private Function f1(x As Single) As Single

f1 = y0 / ((x0 + 1) * Exp(-x0)) * (x + 1) * Exp(-x)

End Function

Private Sub Eiler()

ReDim x(n)

ReDim j(n)

j(0) = y0

For i = 0 To n

x(i) = x0 + h * i

Next i

For i = 0 To n - 1

j(i + 1) = j(i) + h * f(x(i), j(i))

MSFlexGrid1.TextMatrix(1, 0) = Str(x0)

MSFlexGrid1.TextMatrix(i + 2, 0) = Str(x(i + 1))

MSFlexGrid1.TextMatrix(i + 2, 1) = Str(j(i + 1))

MSFlexGrid1.TextMatrix(1, 1) = Str(j(0))

Next i

End Sub

Private Sub Runge()

ReDim y(n)

y(0) = y0

For i = 0 To n

x(i) = x0 + h * i

Next i

For i = 0 To n - 1

k1 = h * f(x(i), y(i))

k2 = h * f(x(i) + h / 2, y(i) + k1 / 2)

k3 = h * f(x(i) + h / 2, y(i) + k2 / 2)

k4 = h * f(x(i) + h, y(i) + k3)

k = (k1 + 2 * k2 + 2 * k3 + k4) / 6

y(i + 1) = y(i) + k

MSFlexGrid1.TextMatrix(1, 2) = Str(y0)

MSFlexGrid1.TextMatrix(i + 2, 2) = Str(y(i + 1))

Next i

End Sub

Private Sub Obhee()

ReDim o(n)

For i = 0 To n

o(0) = y0

x(i) = x0 + h * i

o(i) = f1(x(i))

MSFlexGrid1.TextMatrix(i + 1, 3) = Str(o(i))

Next i

End Sub

Private Sub Command1_Click()

x0 = Val(Text1.Text)

xk = Val(Text2.Text)

h = Val(Text4.Text)

y0 = Val(Text3.Text)

n = (xk - x0) / h

Label6.Caption = Str(x0)

Label5.Caption = Str(xk)

MSFlexGrid1.Rows = n + 2

MSFlexGrid1.TextMatrix(0, 0) = "x"

MSFlexGrid1.TextMatrix(0, 1) = "Ýéëåð"

MSFlexGrid1.TextMatrix(0, 2) = "Ðóíãå-Êóòò"

MSFlexGrid1.TextMatrix(0, 3) = "Îáùåå ðåøåíèå"

Eiler

Runge

Obhee

max = y0

min = y0

For i = 0 To n

If j(i) > max Then

max = j(i)

End If

If j(i) < min Then

min = j(i)

End If

If y(i) > max Then

max = y(i)

End If

If y(i) < min Then

min = y(i)

End If

If o(i) > max Then

max = o(i)

End If

If o(i) < min Then

min = o(i)

End If

Next i

Label4.Caption = Str(max)

Label7.Caption = Str(min)

kx = (6600 - 720) / (xk - x0)

ky = (6600 - 1120) / (max - min)

Picture1.Cls

For i = 1 To n - 1

X1 = 720 + Round(kx * (x(i - 1) - x0))

X2 = 720 + Round(kx * (x(i) - x0))

Y1 = 6600 - Round(ky * (j(i - 1) - min))

Y2 = 6600 - Round(ky * (j(i) - min))

Picture1.Line (X1, Y1)-(X2, Y2), RGB(0, 200, 0)

X1 = 720 + Round(kx * (x(i - 1) - x0))

X2 = 720 + Round(kx * (x(i) - x0))

Y1 = 6600 - Round(ky * (y(i - 1) - min))

Y2 = 6600 - Round(ky * (y(i) - min))

Picture1.Line (X1, Y1)-(X2, Y2), RGB(500, 70, 90)

X1 = 720 + Round(kx * (x(i - 1) - x0))

X2 = 720 + Round(kx * (x(i) - x0))

Y1 = 6600 - Round(ky * (o(i - 1) - min))

Y2 = 6600 - Round(ky * (o(i) - min))

Picture1.Line (X1, Y1)-(X2, Y2), RGB(400, 100, 12)

Next i

End Sub

Private Sub Command2_Click()

End

End Sub

Private Sub Text3_Change()

End Sub

11111

  1.  Решение задачи в MathCad

Эйлер

Рунге-Кутт


Общее

Заключение

 Из двух методов (Эйлера и Рунге-Кутта) по полученным результатам точнее (сравнивая с общим решением) оказался метод Рунге-Кутта. Это объясняется тем что, ведь в отличие от метода Эйлера в методе Рунге-Кутта шаг делится не на 4 отрезка, в результате чего погрешность метода становится меньше.  

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


y

A

e

B

y=y(x)

yi+1

h

yi

x

O

хi

xi+1

α

Eiler(X0, Xk, Y0, N, Y)

h = (Xk – X0)/N

i = 0, …, N - 1

x = X0 + i ∙ h

Yi+1 = Yi + h ∙ F(x, Yi)

End

Rynge4(X0, Xk, Y0, N, Y)

h = (Xk – X0)/N

i = 0, …, N-1

x = X0 + i ∙ h

            k1=h*F(x,Yi)

k2=h*F(x+h/2, Yi +k1/2)

k3= h*F(x+h/2, Yi +k2/2)

k4= h*F(x+h, Yi +k3)

=

k=(k1+2*k2+2*k3+k4)/6

Yi+1= Yi+k

End

      Eiler

ReDim x(n)

ReDim j(n)

j(0) = y0

i = 0, …, n

x(i) = x0 + h * i

i = 0, …, n-1

j(i + 1) = j(i) + h * f(x(i), j(i))

MSFlexGrid1.TextMatrix(1, 0) = Str(x0)

MSFlexGrid1.TextMatrix(i + 2, 0) = Str(x(i + 1))

MSFlexGrid1.TextMatrix(i + 2, 1) = Str(j(i + 1))

MSFlexGrid1.TextMatrix(1, 1) = Str(j(0))

Конец

Runge

ReDim y(n)

g(0) = y0

i = 0, …, n

x(i) = x0 + h * i

i = 0, …, n-1

k1 = h * f(x(i), y(i))

k2 = h * f(x(i) + h / 2, y(i) + k1 / 2)

k3 = h * f(x(i) + h / 2, y(i) + k2 / 2)

k4 = h * f(x(i) + h, y(i) + k3)

k = (k1 + 2 * k2 + 2 * k3 + k4) / 6

y(i + 1) = y(i) + k

MSFlexGrid1.TextMatrix(i + 2, 2) = Str(y(i + 1))

MSFlexGrid1.TextMatrix(1, 2) = Str(y(0))

Конец

Obhee

ReDim o(n)

o(0) = y0

i = 0, …, n

x(i) = x0 + h * i

o(i) = f1(x(i))

Конец

MSFlexGrid1.TextMatrix(i + 1, 3) = Str(o(i))

f1(x)

Конец

f1 = y0 / ((x0 + 1) * Exp(-x0)) * (x + 1) * Exp(-x)

f(a,b)

f = -a * b / (a + 1)

Конец

Начало

y0, x0,xk,h

n = (xk - x0) / h

max = y0

min = y0

MSFlexGrid1.Rows = n + 2

MSFlexGrid1.TextMatrix(0, 0) = "x"

MSFlexGrid1.TextMatrix(0, 1) = "Эйлер "

MSFlexGrid1.TextMatrix(0, 2) = "Рунге-Кутта"

MSFlexGrid1.TextMatrix(0, 3) = "Общее решение"

Label6.Caption = Str(x0)

Label5.Caption = Str(xk)

Eiler

Runge

Obshee

i = 0, … n

j(i) > max

j(i)=max

j(i)=min

y(i) > max

y(i)=max

o(i) > max

o(i)=max

o(i) <min

o(i)=min

y(i)=min

y(i)<min

Label4.Caption = Str(max)

Label7.Caption = Str(min)

kx = (6600 - 720) / (xk - x0)

ky = (6600 - 1120) / (max - min)

= 1, …, n-1

X1 = 720 + Round(kx * (x(i - 1) - x0))

X2 = 720 + Round(kx * (x(i) - x0))

Y1 = 6600 - Round(ky * (j(i - 1) - min))

Y2 = 6600 - Round(ky * (j(i) - min))

X1 = 720 + Round(kx * (x(i - 1) - x0))

X2 = 720 + Round(kx * (x(i) - x0))

Y1 = 6600 - Round(ky * (y(i - 1) - min))

Y2 = 6600 - Round(ky * (y(i) - min))

X1 = 720 + Round(kx * (x(i - 1) - x0))

X2 = 720 + Round(kx * (x(i) - x0))

Y1 = 6600 - Round(ky * (o(i - 1) - min))

Y2 = 6600 - Round(ky * (o(i) - min))

Picture1.Line (X1, Y1)-(X2, Y2), RGB(0, 200, 0)

Picture1.Line (X1, Y1)-(X2, Y2), RGB(500, 70, 90)

Picture1.Line (X1, Y1)-(X2, Y2), RGB(400, 100, 12)

Конец

Command1

Labe31

Text2

Text1

Command2

MSFlexGrid16

Labe21

Text4

Text3

Labe71

Label4

Labe91

Label10

Labe41

Label1

Picture1

Label6

Labe51


 

А также другие работы, которые могут Вас заинтересовать

80670. Коммуникационная политика в маркетинге 76.5 KB
  Планирование и контроль мероприятий ФОССТИС. Цель задачи и правила ФОССТИС Воздействие на покупателя заключающееся в убеждении последнего в том что приобретение именно данного товара выгодно и целесообразно осуществляется различными методами: рекламными посланиями проспекты объявления телевизионные фильмы и т. Все эти средства называются КОММУНИКАЦИОННОЙ ПОЛИТИКОЙ и известны в специальной литературе как средства ФОССТИС формирования спроса и стимулирования сбыта.
80671. Организация сбытовой сети и системы товародвижения в маркетинге 53 KB
  В понятие сбыта включаются факторы: транспортировка складирование хранение доработка продвижение к розничным и оптовым торговым звеньям предпродажная подготовка собственно продажа.Проблемы эффективности рыночного поведения и развития предприятия Наиболее тесно предприятие соприкасается с потребителем в системе сбыта. Руководить предприятием эффективно быть постоянно ориентированным на нужды потребителя СИСТЕМА СБЫТА комплекс состоящий из сбытовой сети предприятия и тех каналов сбыта которые ею пользуются для продажи товаров....
80675. Планирование культуры экономической организации 186.5 KB
  Планирование организационной культуры - вид планирования, в определенном смысле противоположный процессу и результатам финансового планирования. Если финансовое планирование имеет дело с предельно точными, конкретными величинами, то планирование культуры связано с наименее определенным, слабо контролируемым элементом внутрифирменной среды. Основу культуры составляют человек, его поведение