Если у вас возникли сложности с курсовой, контрольной, дипломной, рефератом, отчетом по практике, научно-исследовательской и любой другой работой - мы готовы помочь.
Предоплата всего
Подписываем
по информатике:
Визуализация численных методов.
Решение обыкновенных дифференциальных уравнений.
студент гр.: МЕ-61,
Устинов А.В.
Проверил:
Минина Е.Е.
Екатеринбург
2007 г.
Содержание:
Введение………………………………………………………………….3
1. Постановка задачи…………………………………………………….4
2. Описание методов решения…………………………………………..5
2. 2. Геометрический смысл задачи…………………………………….5
2. 3. Численные методы решения задачи Коши……………………….6
2. 6. Решение поставленной задачи методами Эйлера и Эйлера модифицированного…………………………………………………………….10
2. 6. 1. Метод Эйлера……………………………………………………10
2. 6. 2. Метод Эйлера модифицированный……………………………12
3. Алгоритм решения задачи…………………………………………...14
3. 1. Алгоритмы подпрограмм.………………………………………....14
3. 1. 1. Подпрограмма метода Эйлера………………………………….14
3. 1. 2 Подпрограмма метода Эйлера модифицированного…………..14
3. 1. 3. Подпрограмма общего решения и поиска максимальных значений x и y……………………………………………………………………15
3. 2. Алгоритм функции…………………………………………………15
3. 3. Алгоритм программы………………………………………………16
4. Форма программы…………………………………………………….28
5. Листинг программы…………………………………………………..19
6. Решение задачи в MathCad…………………………………………..21
Заключение………………………………………………………………23
Введение.
В настоящее время существует множество технических систем и технологических процессов, характеристики которых непрерывно меняются со временем. Такие явления обычно подчиняются физическим законам, которые формулируются в виде дифференциальных уравнений.
Дифференциальным называется уравнение, содержащие одну или несколько производных. Лишь очень немногие из них удается решить без помощи вычислительной техники. Поэтому численные методы решения дифференциальных уравнений играют огромную роль в практике инженерных расчетов.
Для решения дифференциальных уравнений удобно использовать языки программирования, так как они позволяют быстро и точно найти решения уравнения и построить график интегральной кривой.
Целью моей работы является решение задачи Коши для дифференциального уравнения 1-го порядка.
Задачами моей работы являются:
1. Изучить численные методы решения дифференциальных уравнений;
2. Самостоятельно вычислить первую точку интегральной кривой заданными методами;
3. Написать программу для построения интегральной кривой;
4. Проверить решение в среде MathCad.
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=2x2+2y |
0 |
1 |
0.1 |
1 |
y=1.5exp(2x)-x2-x-c |
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. Решение поставленной задачи методами Эйлера и Эйлера модифицированного.
2. 6. 1. Метод Эйлера.
1. Строим оси координат;
2. Отмечаем A(1; 2) – первую точку интегральной кривой;
3. Ищем угол наклона касательной к графику в точке A:
α0 = 80°;
4. Строим касательную l0 в точке А под углом α0;
5. Находим х1 по формуле: xi = х0 + ih, где h – шаг интегрирования
x1 = 1 + 1 · 0,05 = 1,05;
6. Проводим прямую x = x1 = 1,05 до пересечения с прямой 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)) = 2 + 0,05 · f(1;2) = 2 + 0,05 · 6 = 2,3
Следовательно, точка B имеет координаты (1,05; 2,3).
Рисунок 8. Решение задачи методом Эйлера.
2. 6. 2. Метод Эйлера модифицированный.
1. Строим оси координат;
2. Отмечаем А(1; 2) – первую точку интегральной кривой;
3. Ищем угол наклона касательной к графику в точке A:
α0 = 80°;
4. Строим касательную l0 в точке А под углом α0;
5. Находим х1 по формуле: xi = х0 + ih, где h – шаг интегрирования
x1 = 1 + 1 · 0,05 = 1,05;
6. Отмечаем середину отрезка x0x1: x0 + h/2, проводим прямую из этой точки до прямой l0, отмечаем точку B(xB; yB);
7. Ищем координаты В:
xB = x0 + h/2 = 1 + 0,05/2 = 1,025
yB = y0 + h/2 · f(x0; y0) = 2 + 0,05/2 · 6 = 2,15
Следовательно, точка B имеет координаты (1,025; 2,15);
8. Ищем угол наклона касательной к графику в точке B:
αB = arctg(f(xB; yB)) = arctg((3 · 2,15)/1,025)) = arctg(6,29) = 1,413рад
αB = 81°;
9. Строим касательную l1 в точке B под углом αB;
10. Проводим прямую x = x1 = 1,05 до пересечения с прямой l1, отмечаем точку C(x1; y1);
11. Ищем y точки C:
y1 = yB + h/2(f(xB;yB)) = 2,15 + 0,05/2 · 6,29 = 2,307
Следовательно, точка C имеет координаты (1,05; 2,307).
Рисунок 9. Решение задачи методом Эйлера модифицированного.
3. Алгоритм решения задачи.
3. 1. Алгоритмы подпрограмм.
3. 1. 1. Подпрограмма метода Эйлера.
3. 1. 2 Подпрограмма метода Эйлера модифицированного.
3. 1. 3. Подпрограмма общего решения и поиска максимальных значений x и y.
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 = (2 * x * x) - 2
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 EilerM()
ReDim x(n + 1)
ReDim em(n + 1)
em(0) = y0
For i = 0 To n
x(i) = Round(x0 + i * h, 3)
em(i + 1) = Round(em(i) + h * f(x(i) + h / 2, em(i) + h / 2 * f(x(i), em(i))), 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)
o(i) = Round(2 * (x(i) ^ 3), 3)
If o(i) > maxy Then maxy = o(i)
If o(i) < miny Then miny = o(i)
If x(i) > maxx Then maxx = x(i)
If x(i) < minx Then minx = x(i)
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) = "y îáù"
MSFlexGrid1.TextMatrix(0, 2) = "y ýéë"
MSFlexGrid1.TextMatrix(0, 3) = "y ýéë ìîä"
Eiler
EilerM
Obhee
For i = 0 To n
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))
Next i
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)
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)
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)
Next i
End Sub
6. Решение задачи в MathCad.
Заключение.
В своей работе я решил дифференциальное уравнение двумя методами: методом Эйлера и методом Эйлера модифицированного. Я произвел вычисления первой точки самостоятельно и построил интегральные кривые в языке программирования Visual Basic и среде Mathcad. Исследуя полученные графики, я увидел, что метод Эйлера модифицированного является более точным методом решения дифференциальных уравнений, так как он дает меньшую погрешность и почти совпадает с графиком общего решения. Также метод Эйлера модифицированного является немного сложнее метода Эйлера.
tg(α) = f(x,y)
α
Эйлер(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
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)
№ докум.
Лист
Изм.
Лист
Дата
Подпись
№ докум.
Лист
Изм.
ε
yi+1
yi
h
α
B
A
0
x
xi+1
xi
y=y(x)
y
C
Δy
α1
α
ε
ε1
xi+1
xi
h
h/2
В
С
А
0
y=y(x)
x
y
Δx
1
1,05
2
2,1
2,2
2,3
l0
Δy
Δx
α0
A
B
1,1
x
y
y
x
1,1
С
A
α0
Δx
Δy
l0
2,3
2,2
2,1
2
1,05
1
B
l1
αB
Начало
y0, x0,xk,h
n = Round((xk - x0) / h)
MSFlexGrid1.Cols = 4
MSFlexGrid1.Rows = n + 2
MSFlexGrid1.TextMatrix(0, 0) = "x"
MSFlexGrid1.TextMatrix(0, 1) = "y общ"
MSFlexGrid1.TextMatrix(0, 2) = "y эйл"
MSFlexGrid1.TextMatrix(0, 3) = "y эйл Эмод"
iler
EilerM
Obhee
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)
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)
i = 1, …, n-1
Picture1.Cls
kx = (Picture1.Width - 1200) / (xk - x0)
ky = (Picture1.Height - 1000) / (maxy - miny)
miny
minx
maxy
maxx
i = 1, …, n-1
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
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
Конец
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
Конец
Конец
x(i) = Round(x0 + i * h, 3)
em(i + 1) = Round(em(i) + h * f(x(i) + h / 2, em(i) + h / 2 * f(x(i), em(i))), 3)
i = 1, …, n
ReDim x(n + 1)
ReDim em(n + 1)
em(0) = y0
EilerM
Конец
x(i) = Round(x0 + (i * h), 3)
o(i) = Round(2 * (x(i) ^ 3), 3)
If o(i) > maxy Then maxy = o(i)
If o(i) < miny Then miny = o(i)
If x(i) > maxx Then maxx = x(i)
If x(i) < minx Then minx = x(i)
i = 1, …, n
ReDim x(n + 1)
ReDim o(n + 1)
maxy = y0
miny = y0
maxx = x0
minx = x0
Начало
Конец
Label1
Labe21
f = (2*x*x)-2
f(a,b)
Labe31
Labe41
Text1
Text2
Text3
Text4
Command1
Label5
Label6
MSFlexGrid16
Picture1
Labe71
Label11
Label10
Labe81
Label12
Labe91
Лист
Дата
Подпись
№ докум.
Лист
Изм.
Лист
Дата
Подпись
№ докум.
Лист
Изм.
Лист
Дата
Подпись
№ докум.
Лист
Изм.
Лист
Дата
Подпись
№ докум.
Лист
Изм.
Лист
Дата
Подпись
№ докум.
Лист
Изм.
Лист
Дата
Подпись
№ докум.
Лист
Изм.
Лист
Дата
Подпись
№ докум.
Лист
Изм.
Лист
Дата
Подпись
№ докум.
Лист
Изм.
Лист
Дата
Подпись
№ докум.
Лист
Изм.
Лист
Дата
Подпись
№ докум.
Лист
Изм.
Лист
Дата
Подпись
№ докум.
Лист
Изм.
Лист
Дата
Подпись
№ докум.
Лист
Изм.
Лист
Дата
Подпись
№ докум.
Лист
Изм.
Лист
Дата
Подпись
№ докум.
Лист
Изм.
Лист
Дата
Подпись
№ докум.
Лист
Изм.
Лист
Дата
Подпись
№ докум.
Лист
Изм.
Лист
Дата
Подпись
№ докум.
Лист
Изм.
Лист
Дата
Подпись
№ докум.
Лист
Изм.
Лист
Дата
Подпись
№ докум.
Лист
Изм.
Лист
Дата
Подпись
№ докум.
Лист
Изм.