Если у вас возникли сложности с курсовой, контрольной, дипломной, рефератом, отчетом по практике, научно-исследовательской и любой другой работой - мы готовы помочь.
Предоплата всего
Подписываем
по информатике:
Визуализация численных методов.
Решение обыкновенных дифференциальных уравнений.
студент гр.:
Проверил:
Минина Е.Е.
Екатеринбург
2007 г.
Содержание:
Введение………………………………………………………………….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
Введение.
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’ + y/x = 3/x |
1 |
1,8 |
0,1 |
0 |
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. 7. Решение поставленной задачи методами Эйлера и Эйлера модифицированного.
2. 7. 1. Метод Эйлера.
1. Строим оси координат;
2. Отмечаем A(1; 2) – первую точку интегральной кривой;
3. Ищем угол наклона касательной к графику в точке A:
4. Строим касательную l0 в точке А под углом α0;
5. Находим х1 по формуле: xi = х0 + ih, где h – шаг интегрирования
x1 = 1 + 1 · 0,1 = 1,1;
6. Проводим прямую x = x1 = 1,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(1;0) = 0 + 0,1 · 3 = 0,3
Следовательно, точка B имеет координаты (1,1; 0,3).
Рисунок 8. Решение задачи методом Эйлера.
2. 7. 2. Метод Эйлера модифицированный.
1. Строим оси координат;
2. Отмечаем А(1; 0) – первую точку интегральной кривой;
3. Ищем угол наклона касательной к графику в точке A:
4. Строим касательную l0 в точке А под углом α0;
5. Находим х1 по формуле: xi = х0 + ih, где h – шаг интегрирования
x1 = 1 + 1 · 0,1 = 1,1;
6. Отмечаем середину отрезка x0x1: x0 + h/2, проводим прямую из этой точки до прямой l0, отмечаем точку B(xB; yB);
7. Ищем координаты В:
xB = x0 + h/2 = 1 + 0,1/2 = 1,05
yB = y0 + h/2 · f(x0; y0) = 0 + 0,1/2 · 3 = 0,15
Следовательно, точка B имеет координаты (1,05; 0,15);
8. Ищем угол наклона касательной к графику в точке B:
αB = arctg(f(xB; yB)) = arctg((3 – 0,15)/1,05)) = arctg(2,71) = 70°
9. Строим касательную l1 в точке B под углом αB;
10. Проводим прямую x = x1 = 1,1 до пересечения с прямой l1, отмечаем точку C(x1; y1);
11. Ищем y точки C:
y1 = yB + h/2(f(xB;yB)) = 0,15 + 0,1/2 · 2,71 = 0,29
Следовательно, точка C имеет координаты (1,1; 0,29).
Рисунок 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 = (3 - b) / a
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(3 * (x(i) - 1) / x(i), 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
EilerM
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)
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.
Заключение.
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)
0,33
0.22
0,11
y
x
A
α0
0
l0
B
α1
α
ε
ε1
xi+1
xi
h
h/2
В
С
А
0
y=y(x)
x
y
y
x
Δx
Δy
A
α0
Δx
Δy
l0
0,33
0.22
0,11
0
1,1
1
1,1
С
1
B
α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 эйл Эмод"
Eiler
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
Obhee
Конец
MSFlexGrid16
Picture1
f = (3 * b) / a
f(a,b)
Эйлер
Эйлер модифицированный
Общее решение
Labe71
Text2
Text1
Labe41
Labe31
Label1
Text3
Labe21
Label11
Text4
Command1
Label10
Labe81
Labe91
Label12
α
xi+1
хi
O
x
yi
h
yi+1
y=y(x)
B
e
A
y