Если у вас возникли сложности с курсовой, контрольной, дипломной, рефератом, отчетом по практике, научно-исследовательской и любой другой работой - мы готовы помочь.
Предоплата всего
Подписываем
Федеральное агентство связи
ГОУ ВПО “СибГУТИ”
Уральский технический институт связи и информатики
по информатике:
Визуализация численных методов.
Решение обыкновенных дифференциальных уравнений.
студентка гр.: СЕ-81
Митрофанова Д.И.
Проверила:
Тюпина О.М.
Екатеринбург
2009 г.
Содержание:
Введение………………………………………………………………….3
1. Постановка задачи…………………………………………………….4
2. Описание методов решения…………………………………………..5
2. 2. Геометрический смысл задачи…………………………………….5
2. 3. Численные методы решения задачи Коши……………………….6
2. 7. Решение поставленной задачи методами Эйлера и Рунге-Кутта…………………………………………………………….12
2. 7. 1. Метод Эйлера……………………………………………………12
2. 7. 2. Метод Рунге-Кутта……………………………13
3. Алгоритм решения задачи…………………………………………...16
3. 1. Алгоритмы подпрограмм.………………………………………....16
3. 1. 1. Подпрограмма метода Эйлера………………………………….16
3. 1. 2 Подпрограмма метода Эйлера модифицированного…………..16
3. 1. 3. Подпрограмма общего решения и поиска максимальных значений x и y……………………………………………………………………17
3. 2. Алгоритм функции…………………………………………………17
3. 3. Алгоритм программы………………………………………………19
4. Форма программы…………………………………………………….20
5. Листинг программы…………………………………………………..21
6. Решение задачи в MathCad…………………………………………..23
Заключение………………………………………………………………25
Введение
При решении различных задач математики, физики, химии и других наук часто пользуются математическими моделями в виде дифференциальных уравнений связывающих независимую переменную, искомую функцию и ее производные. Например, исследуя полученные дифференциальные уравнения вместе с дополнительными условиями, мы можем получить сведения о происходящем явлении, иногда может узнать его прошлое и будущее. Изучение математической модели математическими методами позволяет не только получить качественные характеристики физических явлений и рассчитать с заданной степенью точности ход реального процесса, но и дает возможность проникнуть в суть физических явлений, а иногда предсказать и новые физические эффекты. Бывает, что сама природа физического явления подсказывает и подходы, и методы математического исследования. Критерием правильности выбора математической модели является практика, сопоставление данных математического исследования с экспериментальными данными. Не всегда удается решить дифференциальное уравнение без помощи компьютера, поэтому создание все лучших и удобных программ для решения уравнений всегда будет являться актуальным вопросом.
Целью данной курсовой работы является решение дифференциального уравнения двумя численными методами: методом Эйлера и методом Рунге-Кутта 4 порядка точности.
Для достижения цели я поставил перед собой следующие задачи:
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+1)=y+2 |
0 |
0,8 |
0,1 |
0 |
y=(x+1)*c-2 |
2. Описание методов решения
Чтобы решить обыкновенное дифференциальное уравнение, необходимо знать значения зависимой переменной и (или) её производных при некоторых значениях независимой переменной. Если эти дополнительные условия задаются при одном значении независимой переменной, то такая задача называется задачей с начальными условиями, или задачей Коши. Часто в задаче Коши в роли независимой переменной выступает время.
Задачу Коши можно сформулировать следующим образом:
Пусть дано дифференциальное уравнение и начальное условие y(x0) = у0. Требуется найти функцию у(x), удовлетворяющую как указанному уравнению, так и начальному условию.
Численное решение задачи Коши сводится к табулированию искомой функции.
График решения дифференциального уравнения называется интегральной кривой.
2. 2. Геометрический смысл задачи.
y’ = f(x,y) - тангенс угла наклона касательной к графику решения в точке (х, у) к оси 0Х, - угловой коэффициент (рис. 1).
Рисунок 1. Геометрический смысл задачи Коши.
Существование решения:
Если правая часть f(x, y) непрерывна в некоторой области R, определяемой неравенствами
|x-x0| < а; |y-y0| < b,
то существует, по меньшей мере, одно решение у = у(х), определённое в окрестности |х – х0| < h, где h - положительное число.
Это решение единственно, если в R выполнено условие Липшица
|f(x,y)-f(x,y)| ≤N|y-y|(x,y),
где N - некоторая постоянная (константа Липшица), зависящая, в общем случае, от а и b. Если f(x, у) имеет ограниченную производную
f’y(x, y) в R, то можно положить N = мах |f’y(х, у)| при (х, y) принадлежащим R.
2. 3. Численные методы решения задачи Коши
При использовании численных методов выполняется замена отрезка [х0, X] - области непрерывного изменения аргумента х множеством . состоящего из конечного числа точек х0 < х1 < ... < xn = Х - сеткой.
При этом xi называют узлами сетки.
Во многих методах используются равномерные сетки с шагом:
Задача Коши, определённая ранее на непрерывном отрезке [х0, X], заменяется её дискретным аналогом - системой уравнений, решая которую можно последовательно найти значения y1, y2,…,yn - приближённые значения функции в узлах сетки.
Численное решение задачи Коши широко применяется в различных областях науки и техники, и число разработанных для него методов достаточно велико. Эти методы могут быть разделены на следующие группы.
Одношаговые методы, в которых для нахождения следующей точки на
кривой у = f(x) требуется информация лишь об одном предыдущем шаге.
Одношаговыми являются метод Эйлера и методы Рунге - Кутта.
Методы прогноза и коррекции (многошаговые), в которых для отыскания следующей точки кривой у = f(x) требуется информация более чем об одной из предыдущих точек. Чтобы получить достаточно точное численное значение, часто прибегают к итерации. К числу таких методов относятся методы Милна, Адамса - Башфорта и Хемминга.
Явные методы, в которых функция Ф не зависит от yn+1.
Неявные методы, в которых функция Ф зависит от yn+1.
Иногда этот метод называют методом Рунге-Кутта первого порядка точности.
Данный метод одношаговый. Табулирование функции происходит поочередно в каждой точке. Для расчета значения функции в очередном узле необходимо использовать значение функции в одном предыдущем узле.
Пусть дано дифференциальное уравнение первого порядка:
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) и . Получим
Отсюда Δу = h ∙ f(xi, yi).
Подставим в это выражение формулу yi+1 = yi + Δy, а затем преобразуем его. В результате получаем формулу расчета очередной точки интегральной функции:
.
Из формулы видно, что для расчета каждой следующей точки интегральной функции необходимо знать значение только одной предыдущей точки. Таким образом, зная начальные условия, можно построить интегральную кривую на заданном промежутке.
Блок-схема процедуры решения дифференциального уравнения методом Эйлера приведена на рисунке 3.
F(x, у) - заданная функция – должна
быть описана отдельно.
Входные параметры:
Х0, XK—начальное и конечное
значения независимой переменной;
Y0 – значение y0 из начального условия
y(x0) = y0;
N - количество отрезков разбиения;
Выходные параметры:
У - массив значений искомого решения
в узлах сетки.
Рисунок 3. Блок-схема процедуры решения дифференциального уравнения методом Эйлера.
Метод Эйлера - один из простейших методов численного решения обыкновенных дифференциальных уравнений. Но существенным его недостатком является большая погрешность вычислений. На рисунке 2 погрешность вычислений для i-гo шага обозначена ε. С каждым шагом погрешность вычислений увеличивается.
Для уменьшения погрешности вычислений часто используется модифицированный метод Эйлера. Этот метод имеет так же следующие названия: метод Эйлера-Коши или метод Рунге-Кутта второго порядка точности.
Пусть дано дифференциальное уравнение первого порядка
с начальным условием:
Выберем шаг h и введём обозначения:
xi = x0 + ih и yi = y(xi), где i = 0, 1, 2, ...,
xi - узлы сетки,
yi - значение интегральной функции в узлах.
При использовании модифицированного метода Эйлера шаг h делится на два отрезка.
Иллюстрации к решению приведены на рисунке 4.
Рисунок 4. Метод Эйлера модифицированный
Проведем решение в несколько этапов:
yi+1 = yi + h ∙ f(xi + h/2, yi + h/2 ∙ f(xi, yi)).
Модифицированный метод Эйлера дает меньшую погрешность. На рисунке 4 это хорошо видно. Так величина εl характеризует погрешность метода Эйлера, а ε - погрешность метода Эйлера модифицированного.
Блок-схема процедуры решения дифференциального уравнения методом Эйлера модифицированным приведена на рисунке 5.
F(x, у) - заданная функция - должна
быть описана отдельно.
Входные параметры:
Х0, XК - начальное и конечное
значения независимой
переменной;
Y0 – значение y0 из начального условия
y(x0)=y0;
N - количество отрезков разбиения;
Выходные параметры:
Y - массив значений искомого решения
в узлах сетки.
Рисунок 5. Блок-схема процедуры решения дифференциального уравнения методом Эйлера модифицированным.
2.6 Метод Рунге-Кутта 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. 7. Решение поставленной задачи методами Эйлера и Рунге-Кутта
2. 7. 1. Метод Эйлера
1. Строим оси координат;
2. Отмечаем A(0; 0) – первую точку интегральной кривой;
3. Ищем угол наклона касательной к графику в точке A:
4. Строим касательную l0 в точке А под углом α0;
5. Находим х1 по формуле: xi = х0 + ih, где h – шаг интегрирования
x1 = 0 + 1 · 0,1 = 0,1;
6. Проводим прямую x = x1 = 0,1 до пересечения с прямой l0, отмечаем точку B(x1; y1);
7. Ищем y точки B:
Из прямоугольного треугольника ABC ,
Δy = y1 – y0,
Δx = x1 – x0 = h,
f(x0; y0) = (y1 – y0)/h =>
y1 = y0 + h · (f(x0; y0)) = 0 + 0,1 · f(0;0) = 0 + 0,1 · 2 = 0,2
Следовательно, точка B имеет координаты (0,1; 0,2).
2.7.2. Метод Рунге-Кутта 4 порядка
1. Строим оси координат;
2. Отмечаем А(0; 0) – первую точку интегральной кривой;
3. Ищем угол наклона касательной к графику в точке A:
4. Строим касательную l0 в точке А под углом α0;
5. Находим х1 по формуле: xi = х0 + ih
x1 = 0 + 1 · 0,1 = 0,1;
k1=0,1·f(0,0)=0,1·2=0,2
k2=0,1· f(0,1+0,1/2;0+0,2/2)= 0,2
k3=0,1· f(0,1+0,1/2;0+0,2/2)= 0,2
k4=0,1· f(0,1+0,1;0+0,2)= 0,2
∆y1=(0,2+2·0,2+2·0,2+0,2)/6=0,2
∆y2=0+0,2=0,2
Следовательно, следующая точка графика решения имеет координаты (0,1; 0,2)
. Алгоритм решения задачи
3. 1. Алгоритмы подпрограмм
3. 1. 1. Подпрограмма метода Эйлера
3.1.2 Подпрограмма метода Рунге-Кутта 4 порядка
3. 1. 3. Подпрограмма общего решения
3. 2. Алгоритм функции
3. 3. Алгоритм программы
4. Форма программы.
5. Листинг программы.
Dim x(), e(), em(), o() As Single
Private i, n As Integer
Private x0, xk, y0, h, miny, maxy, minx, maxx As Single
Function f(a, b) As Single
f = (b + 2) / (a + 1)
End Function
Private Sub Eiler()
ReDim x(n + 1)
ReDim e(n + 1)
e(0) = y0
For i = 0 To n
x(i) = Round(x0 + (i * h), 3)
e(i + 1) = Round(e(i) + h * f(x(i), e(i)), 3)
Next i
End Sub
Private Sub RungeKutta()
ReDim x(n + 1)
ReDim em(n + 1)
em(0) = y0
For i = 0 To n
x(i) = Round(x0 + i * h, 3)
k1 = h * f(x(i), em(i))
k2 = h * f(x(i) + (h / 2), em(i) + (k1 / 2))
k3 = h * f(x(i) + (h / 2), em(i) + (k2 / 2))
k4 = h * f(x(i) + h, em(i) + k3)
k = (k1 + 2 * k2 + 2 * k3 + k4) / 6
em(i + 1) = Round(em(i) + k, 3)
Next i
End Sub
Private Sub Obhee()
ReDim x(n + 1)
ReDim o(n + 1)
maxy = y0
miny = y0
maxx = x0
minx = x0
For i = 0 To n
x(i) = Round(x0 + i * h, 3)
c = (y0 + 2) / (x0 + 1)
o(i) = Round((x(i) + 1) * c - 2, 3)
Next i
End Sub
Private Sub Command1_Click()
x0 = Val(Text1.Text)
y0 = Val(Text2.Text)
xk = Val(Text3.Text)
h = Val(Text4.Text)
n = Round((xk - x0) / h)
MSFlexGrid1.Cols = 4
MSFlexGrid1.Rows = n + 2
MSFlexGrid1.TextMatrix(0, 0) = "x"
MSFlexGrid1.TextMatrix(0, 1) = "Общее рещение"
MSFlexGrid1.TextMatrix(0, 2) = "эйлер"
MSFlexGrid1.TextMatrix(0, 3) = "Рунге-Кутт"
Eiler
RungeKutta
Obhee
For i = 0 To n
MSFlexGrid1.TextMatrix(i + 1, 0) = Str(x(i))
MSFlexGrid1.TextMatrix(i + 1, 1) = Str(o(i))
MSFlexGrid1.TextMatrix(i + 1, 2) = Str(e(i))
MSFlexGrid1.TextMatrix(i + 1, 3) = Str(em(i))
Next i
minx = x(0)
maxx = x(n)
miny = o(0)
maxy = o(n)
If e(n) > o(n) Then maxy = e(n)
If em(n) > o(n) Then maxy = em(n)
If e(n) > em(n) Then maxy = e(n)
Label10.Caption = Str(miny)
Label11.Caption = Str(maxy)
Label8.Caption = Str(minx)
Label12.Caption = Str(maxx)
Picture1.Cls
kx = (Picture1.Width - 1200) / (xk - x0)
ky = (Picture1.Height - 1000) / (maxy - miny)
For i = 0 To n - 1
z1 = Round(720 + (x(i) - x0) * kx)
z2 = Round(5400 - (e(i) - miny) * ky)
z3 = Round(720 + (x(i + 1) - x0) * kx)
z4 = Round(5400 - (e(i + 1) - miny) * ky)
Picture1.Line (z1, z2)-(z3, z4), RGB(0, 0, 9999)
Next i
For i = 0 To n - 1
z1 = Round(720 + (x(i) - x0) * kx)
z2 = Round(5400 - (em(i) - miny) * ky)
z3 = Round(720 + (x(i + 1) - x0) * kx)
z4 = Round(5400 - (em(i + 1) - miny) * ky)
Picture1.Line (z1, z2)-(z3, z4), RGB(0, 9999, 0)
Next i
For i = 0 To n - 1
z1 = Round(720 + (x(i) - x0) * kx)
z2 = Round(5400 - (o(i) - miny) * ky)
z3 = Round(720 + (x(i + 1) - x0) * kx)
z4 = Round(5400 - (o(i + 1) - miny) * ky)
Picture1.Line (z1, z2)-(z3, z4), RGB(9999, 0, 0)
Next i
End Sub
6. Решение задачи в MathCad.
Метод Эйлера
Метод Эйлера IV порядка точности (Метод Рунге-Кутта)
Общее решение
Заключение
При расчете уравнения ,двумя методами (Эйлера и Рунге-Кутта), получил значения сходные с общим, хотя метод Рунге-Кутта является наиболее точным. Это совпадение обуславливается маленьким шагом и небольшим диапазоном конечных значений.
В ходе работы я выполнил все поставленные задачи: написал программу для решения данного дифференциального уравнения двумя численными методами в программе Visual Basic, проверил решение с помощью приложения MathCAD,сравнил полученные разными методами результаты с общим решением. Поэтому я считаю что полностью выполнил поставленное передо мной задание курсовой работы.
tg(α) = f(x,y)
α
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
End
Yi+1 = Yi + h ∙ F(x + h/2, Yi + h/2 ∙ F(xi, yi))
x = X0 + i ∙ h
i = 0, …, N-1
h = (Xk – X0)/N
EilerM(X0, Xk, Y0, N, Y)
α1
α
ε
ε1
xi+1
xi
h
h/2
В
С
А
0
y=y(x)
x
y
em(i + 1) = Round(em(i) + k, 3)
k2 = h * f(x(i) + (h / 2), em(i) + (k1 / 2))
k3 = h * f(x(i) + (h / 2), em(i) + (k2 / 2))
k4 = h * f(x(i) + h, em(i) + k3)
k = (k1 + 2 * k2 + 2 * k3 + k4) / 6
x(i) = Round(x0 + i * h, 3)
k1 = h * f(x(i), em(i))
h=(xk-x0)/n
Конец
i = 1, …, n
Конец
z1 = Round(720 + (x(i) - x0) * kx)
z2 = Round(5400 - (o(i) - miny) * ky)
z3 = Round(720 + (x(i + 1) - x0) * kx)
z4 = Round(5400 - (o(i + 1) - miny) * ky)
Picture1.Line (z1, z2)-(z3, z4)
i = 1, …, n-1
RungeKutta
z1 = Round(720 + (x(i) - x0) * kx)
z2 = Round(5400 - (em(i) - miny) * ky)
z3 = Round(720 + (x(i + 1) - x0) * kx)
z4 = Round(5400 - (em(i + 1) - miny) * ky)
Picture1.Line (z1, z2)-(z3, z4)
End
Yi+1= Yi+k
k=(k1+2*k2+2*k3+k4)/6
k4= h*F(x+h, Yi +k3)
=
k3= h*F(x+h/2, Yi +k2/2)
k2=h*F(x+h/2, Yi +k1/2)
k1=h*F(x,Yi)
x = X0 + i ∙ h
i = 0, …, N-1
i = 1, …, n-1
z1 = Round(720 + (x(i) - x0) * kx)
z2 = Round(5400 - (e(i) - miny) * ky)
z3 = Round(720 + (x(i + 1) - x0) * kx)
z4 = Round(5400 - (e(i + 1) - miny) * ky)
Picture1.Line (z1, z2)-(z3, z4)
h = (Xk – X0)/N
i = 1, …, n-1
Picture1.Cls
kx = (Picture1.Width - 1200) / (xk - x0)
ky = (Picture1.Height - 1000) / (maxy - miny)
miny
minx
maxy
maxx
MSFlexGrid1.TextMatrix(i + 1, 2) = Str(e(i))
MSFlexGrid1.TextMatrix(i + 1, 3) = Str(em(i))
MSFlexGrid1.TextMatrix(i + 1, 1) = Str(o(i))
MSFlexGrid1.TextMatrix(i + 1, 0) = Str(x(i))
i = 1, …, n
Eiler
RungeKutta
Obhee
= Round((xk - x0) / h)
MSFlexGrid1.Cols = 4
MSFlexGrid1.Rows = n + 2
MSFlexGrid1.TextMatrix(0, 0) = "x"
MSFlexGrid1.TextMatrix(0, 1) = "Эйлер"
MSFlexGrid1.TextMatrix(0, 2) = "Рунге-Кутта"
MSFlexGrid1.TextMatrix(0, 3) = "Общее реш."
y0, x0,xk,h
Начало
Конец
f = (y+2)/(x+1)
f(x,y)
x(i) = Round(x0 + i * h, 3)
c = (y0 + 2) / (x0 + 1)
Eiler
ReDim x(n + 1)
ReDim e(n + 1)
e(0) = y0
x(i) = Round(x0 + (i * h), 3)
e(i + 1) = Round(e(i) + h * f(x(i), e(i)), 3)
i = 1, …, n
Конец
Конец
i = 1, …, n
o(i) = Round((x(i) + 1) * c - 2, 3)
maxy = y0
miny = y0
maxx = x0
minx = x0
Obhee
MSFlexGrid16
Picture1
Labe71
Text2
Text1
Labe41
Labe31
Label1
Text3
Labe21
Label11
Rynge4(X0, Xk, Y0, N, Y)
Text4
Command1
Label10
Labe81
Labe91
Label12
α
xi+1
хi
O
x
yi
h
yi+1
y=y(x)
B
e
A
y