Если у вас возникли сложности с курсовой, контрольной, дипломной, рефератом, отчетом по практике, научно-исследовательской и любой другой работой - мы готовы помочь.
Предоплата всего
Подписываем
Федеральное агентство связи
ГОУ ВПО «Сибирский государственный университет телекоммуникаций и информатики»
Уральский технический институт связи и информатики (филиал)
Визуализация численных методов.
Решение обыкновенных дифференциальных уравнений.
Выполнила: студентка гр. МЕ-71
Долматова. К. А
Проверила: Минина Е.Е.
Содержание:
Введение………………………………………………………………….3
1. Постановка задачи…………………………………………………….4
2. Описание методов решения…………………………………………..5
2. 2. Геометрический смысл задачи…………………………………….5
2. 3. Численные методы решения задачи Коши……………………….6
2. 6. Решение поставленной задачи методами Эйлера и Рунге-Кутта 4го порядка…………………………………………………………….11
2. 6. 1. Метод Эйлера ……….……………………11
2. 6. 2. Метод Рунге-Кутта 4-го порядка ……………………………12
3. Алгоритм решения задачи…………………………………………...13
3. 1. Алгоритмы подпрограмм.………………………………………....13
3. 1. 1. Подпрограмма метода Эйлера ……………………..13
3. 1. 2. Подпрограмма метода Рунге-Кутта 4-го порядка ………..14
3. 1. 3. Подпрограмма общего решения ………..…………………15
3. 2. Алгоритм функции…………………………………………………15
3. 3. Алгоритм программы…………………………………16-18
4. Форма программы………………………………………………19-20
5. Листинг программы………………………………………………21-23
6. Решение задачи в MathCad…………………………………………..24
Заключение………………………………………………………………25
2
Курсовая работа является первой самостоятельной работой в процессе моего обучения. Она ставит своей целью закрепить знания по дисциплине “Информатика”. Тема данной курсовой “Визуализация численных методов. Решение обыкновенных дифференциальных уравнений”.
Дифференциальным называется уравнение, содержащее одну или несколько производных.
Численные методы решения дифференциальных уравнений играют огромную роль в практике инженерных расчетов.
Для решения таких уравнений лучше всего использовать языки программирования, так как они более точно могут найти решение того или иного уравнения.
Целью данной работы является решение задачи Коши для дифференциального уравнения 1-го порядка.
Основные задачи, которые я для себя поставила это:
1.Познакомится с численными методами решения дифференциального уравнения(Метод Эйлера и Рунге-Кутта 4го порядка)
2.На основе данных методов найти точки для построения графика
3. Написать программу на языке программирования VusuaiBasic;
4. Проверить решение задачи в 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=-0.05*y |
1 |
2 |
0.1 |
2 |
y=c*exp(- 0.05*x) |
4
2. Описание методов решения
Чтобы решить обыкновенное дифференциальное уравнение, необходимо знать значения зависимой переменной и (или) её производных при некоторых значениях независимой переменной. Если эти дополнительные условия задаются при одном значении независимой переменной, то такая задача называется задачей с начальными условиями, или задачей Коши. Часто в задаче Коши в роли независимой переменной выступает время.
Задачу Коши можно сформулировать следующим образом:
Пусть дано дифференциальное уравнение и начальное условие y(x0) = у0. Требуется найти функцию у(x), удовлетворяющую как указанному уравнению, так и начальному условию.
Численное решение задачи Коши сводится к табулированию искомой функции.
График решения дифференциального уравнения называется интегральной кривой.
2. 2. Геометрический смысл задачи
y’ = f(x,y) - тангенс угла наклона касательной к графику решения в точке (х, у) к оси 0Х, - угловой коэффициент (рис. 1).
Рисунок 1. Геометрический смысл задачи Коши.
Существование решения:
Если правая часть f(x, y) непрерывна в некоторой области R, определяемой неравенствами
|x-x0| < а; |y-y0| < b,
5
то существует, по меньшей мере, одно решение у = у(х), определённое в окрестности |х – х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.
6
Иногда этот метод называют методом Рунге-Кутта первого порядка точности.
Данный метод одношаговый. Табулирование функции происходит поочередно в каждой точке. Для расчета значения функции в очередном узле необходимо использовать значение функции в одном предыдущем узле.
Пусть дано дифференциальное уравнение первого порядка:
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, а затем преобразуем его. В результате получаем формулу расчета очередной точки интегральной функции:
.
7
Из формулы видно, что для расчета каждой следующей точки интегральной функции необходимо знать значение только одной предыдущей точки. Таким образом, зная начальные условия, можно построить интегральную кривую на заданном промежутке.
Блок-схема процедуры решения дифференциального уравнения методом Эйлера приведена на рисунке 3.
F(x, у) - заданная функция – должна
быть описана отдельно.
Входные параметры:
Х0, XK—начальное и конечное
значения независимой переменной;
Y0 – значение y0 из начального условия
y(x0) = y0;
8
N - количество отрезков разбиения;
Выходные параметры:
У - массив значений искомого решения
в узлах сетки.
Рисунок 3. Блок-схема процедуры решения дифференциального уравнения методом Эйлера.
Метод Эйлера - один из простейших методов численного решения обыкновенных дифференциальных уравнений. Но существенным его недостатком является большая погрешность вычислений. На рисунке 2 погрешность вычислений для i-гo шага обозначена ε. С каждым шагом погрешность вычислений увеличивается.
2.5 Метод Рунге-Кутта 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 порядка точности.
9
Блок-схема процедуры решения дифференциального уравнения методом Рунге-Кутта приведена на рисунке 6.
F(x, у) - заданная функция - должна
быть описана отдельно.
Входные параметры:
Х0, XК - начальное и конечное
значения независимой
переменной;
Y0 – значение y0 из начального условия
y(x0)=y0;
N - количество отрезков разбиения;
Выходные параметры:
Y - массив значений искомого решения
в узлах сетки.
10
2.6 Решение поставленной задачи методами Эйлера и Рунге-Кутта 4го порядка
2.6.1 Расчет первой точки методом Эйлера
1.Строим оси координат
2.Отмечаем начальную точку интегральной кривой А(х0,у0)
3.Ищем угол наклона касательной к графику в точке А(1;2)
y’(x0) = f(x0;y0)
α0 = arctg(f(x0;y0))
α0 = arctg((-0.05*2)) = arctg(-0.1)
4.Строим касательную l0 в т. А под углом
5.Находим х1 по формуле:
x1=x0+i*h
x1=1+1*0.1=1.1
6. Проводим прямую х=хi=0.2 до пересечения отмеченной точки В(xi,yi) c прямой l0
7.Ищем у у т.В
x=xi-x0=h
f(x0,y0)=(yi-y0)*hyi=y0+h*(f(x0,y0))=2+0.1*1=2.1 т.В имеет координаты (1.1,2.1)
11
1. Строим оси координат;
2. Отмечаем А(1; 2) – первую точку интегральной кривой;
3. Ищем угол наклона касательной к графику в точке A:
4. Строим касательную l0 в точке А под углом α0;
5. Находим х1 по формуле:
x1=x0+i*h
x1=1+1*0.1=1.1
Находим по формулам:
K1=0.1*(1,2)=-0.05*2=-0.1
K2=0.1*((1+0.1)/2;2+(-0.1))=1.05*1.9=0.1995
K3=0.1*(1+0.1/2;2+0.1995)=0.2301
K4=0.1(1+0.1/2;2+0.2301)=0.2341
y=k=(-0.1+2*0.1995+2*0.2301+0.2341)/6=0.242
т.В(0.1;0,242)
12
3. Алгоритм решения задачи
3. 1. Алгоритмы подпрограмм
3. 1. 1 Подпрограмма метода Эйлера
13
3.1.2 Подпрограмма метода Рунге-Кутта 4 порядка
14
3. 1. 3. Подпрограмма общего решения
3. 2. Алгоритм функции
15
3. 3. Алгоритм программы
16
17
18
4. Форма программы
19
20
5. Листинг программы
Dim x() As Single
Dim y() As Single
Dim g() As Single
Dim t() 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 = -0.05 * b * a
End Function
Private Function v(x As Single) As Single
v = (y0 / Exp(-0.05 * x0)) * Exp(-0.05 * x)
End Function
Private Sub Eiler()
ReDim x(n)
ReDim y(n)
y(0) = y0
For i = 0 To n
x(i) = x0 + h * i
Next i
For i = 0 To n - 1
y(i + 1) = y(i) + h * f(x(i), y(i))
MSFlexGrid1.TextMatrix(1, 0) = Str(x0)
MSFlexGrid1.TextMatrix(i + 2, 0) = Str(x(i + 1))
MSFlexGrid1.TextMatrix(1, 1) = Str(y0)
MSFlexGrid1.TextMatrix(i + 2, 1) = Str(y(i + 1))
Next i
End Sub
Private Sub Rynge_Kytt()
ReDim g(n)
g(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), g(i))
k2 = h * f(x(i) + h / 2, g(i) + k1 / 2)
k3 = h * f(x(i) + h / 2, g(i) + k2 / 2)
k4 = h * f(x(i) + h, g(i) + k3)
k = (k1 + 2 * k2 + 2 * k3 + k4) / 6
g(i + 1) = g(i) + k
MSFlexGrid1.TextMatrix(i + 2, 2) = Str(g(i + 1))
MSFlexGrid1.TextMatrix(1, 2) = Str(g(0))
Next i
End Sub
Private Sub Obshee()
ReDim t(n)
For i = 0 To n
t(0) = y0
x(i) = x0 + h * i
t(i) = v(x(i))
MSFlexGrid1.TextMatrix(i + 1, 3) = Str(t(i))
Next i
End Sub
Private Sub Command1_Click()
x0 = Val(Text1.Text)
xk = Val(Text2.Text)
h = Val(Text3.Text)
y0 = Val(Text4.Text)
n = (xk - x0) / h
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
Rynge_Kytt
Obshee
max = y0
min = y0
For i = 0 To n
If y(i) > max Then
max = y(i)
End If
If y(i) < min Then
min = y(i)
End If
If g(i) > max Then
max = g(i)
End If
22
If g(i) < min Then
min = g(i)
End If
If t(i) > max Then
max = t(i)
End If
If t(i) < min Then
min = t(i)
End If
Next i
Label4.Caption = Str(max)
Label7.Caption = Str(min)
kx = (6600 - 960) / (xk - x0)
ky = (5160 - 240) / (max - min)
Picture1.Cls
For i = 1 To n - 1
X1 = 960 + Round(kx * (x(i - 1) - x0))
X2 = 960 + Round(kx * (x(i) - x0))
Y1 = 5160 - Round(ky * (y(i - 1) - min))
Y2 = 5160 - Round(ky * (y(i) - min))
Picture1.Line (X1, Y1)-(X2, Y2), RGB(9999, 0, 0)
X1 = 960 + Round(kx * (x(i - 1) - x0))
X2 = 960 + Round(kx * (x(i) - x0))
Y1 = 5160 - Round(ky * (g(i - 1) - min))
Y2 = 5160 - Round(ky * (g(i) - min))
Picture1.Line (X1, Y1)-(X2, Y2), RGB(0, 9999, 0)
X1 = 960 + Round(kx * (x(i - 1) - x0))
X2 = 960 + Round(kx * (x(i) - x0))
Y1 = 5160 - Round(ky * (t(i - 1) - min))
Y2 = 5160 - Round(ky * (t(i) - min))
Picture1.Line (X1, Y1)-(X2, Y2), RGB(0, 0, 9999)
Next i
End Sub
Private Sub Command2_Click()
End
End Sub
Private Sub Form_Load()
End Sub
23
6. Решение задачи в MathCad
24
Заключение
По написанию данной курсовой работы можно сделать некоторые выводы. В ней я описала методы Эйлера и Рунге-Кутта четвертого порядка, с помощью которых была решена задача Коши первого порядка.
Были задействованы такие программы, как Visual Basic и Math Cad, использовались построения блок-схем основных процедур, разработка программного кода(Листинг программы).
Также мы закрепили знания по изученным нам уже ранее программам.
В принципе, данная работа была очень интересной ,но относящаяся к среднему уровню сложности для меня, так как безусловно, требующая большого количества времени и труда.
25
α
tg(α) = f(x,y)
x = X0 + i ∙ h
i = 0, …, N - 1
h = (Xk – X0)/N
Eiler(X0, Xk, Y0, N, Y)
α
xi+1
хi
O
x
yi
h
yi+1
y=y(x)
B
e
A
y
Yi+1 = Yi + h ∙ F(x, Yi)
End
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
h = (Xk – X0)/N
Rynge4(X0, Xk, Y0, N, Y)
Конец
MSFlexGrid1.TextMatrix(1, 0) = Str(x0)
MSFlexGrid1.TextMatrix(i + 2, 0) = Str(x(i + 1))
MSFlexGrid1.TextMatrix(1, 1) = Str(y0)
MSFlexGrid1.TextMatrix(i + 2, 1) = Str(y(i + 1))
Конец
y(i + 1) = y(i) + h * f(x(i),y(i))
i = 0, …, n-1
x(i) = x0 + h * i
i = 0, …, n
ReDim x(n )
ReDim y(n )
y(0) = y0
Eiler
Конец
MSFlexGrid1.TextMatrix(i + 2, 2) = Str(g(i + 1))
MSFlexGrid1.TextMatrix(1, 2) = Str(g(0))
k1 = h * f(x(i), g(i))
k2 = h * f(x(i) + h / 2, g(i) + k1 / 2)
k3 = h * f(x(i) + h / 2, g(i) + k2 / 2)
k4 = h * f(x(i) + h, g(i) + k3)
k = (k1 + 2 * k2 + 2 * k3 + k4) / 6
g(i + 1) = g(i) + k
i = 0, …, n-1
x(i) = x0 + h * i
i = 0, …, n
ReDim g(n)
g(0) = y0
Rynge_Kytt
MSFlexGrid1.TextMatrix(i + 1, 3) = Str(t(i))
Конец
x(i) = x0 + h * i
t(i) = v(x(i))
i = 0, …, n
ReDim t(n)
t(0) = y0
Obshee
Конец
v = v = (y0 / Exp(-0.05 * x0)) * Exp(-0.05 * x)
V(x)
Конец
f = -0.05*y
f(x,y)
y(i)<min
g(i)=max
g(i) > max
y(i)=min
y(i)=max
y(i) > max
i = 0, … n
Eiler
Rynge_Kytt
Obshee
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)
n = (xk - x0) / h
max = y0
min = y0
y0, x0,xk,h
Начало
X1 = 960 + Round(kx * (x(i - 1) - x0))
X2 = 960 + Round(kx * (x(i) - x0))
Y1 = 5160 - Round(ky * (t(i - 1) - min))
Y2 = 5160 - Round(ky * (t(i) - min))
X1 = 960 + Round(kx * (x(i - 1) - x0))
X2 = 960 + Round(kx * (x(i) - x0))
Y1 = 5160 - Round(ky * (g(i - 1) - min))
Y2 = 5160 - Round(ky * (g(i) - min))
X1 = 960 + Round(kx * (x(i - 1) - x0))
X2 = 960 + Round(kx * (x(i) - x0))
Y1 = 5160 - Round(ky * (y(i - 1) - min))
Y2 = 5160 - Round(ky * (y(i) - min))
i = 1, …, n-1
kx = (6600 - 960) / (xk - x0)
ky = (5160 - 240) / (max - min)
Label4.Caption = Str(max)
Label7.Caption = Str(min)
g(i)<min
(i)=min
t(i)=min
t(i) <min
t(i)=max
t(i) > max
Конец
Picture1.Line (X1, Y1)-(X2, Y2), RGB(100, 200, 200)
Picture1.Line (X1, Y1)-(X2, Y2), RGB(400, 300, 300)
Picture1.Line (X1, Y1)-(X2, Y2), RGB(500, 100, 100)
MSFlexGrid1
Picture 1
Label10
Label5
Label6
Label7
Label9
Label4
Command1
Command2
Label1
Label2
Label3
Label8
Text1
Text2
Text3
Text4