Реализация и примеры
В рабочей книге couchy.xls представлен вариант реализации рассмотренных методов решения задачи Коши.
Задание правой части дифференциального уравнения y'=f(x, y), начальной точки, интервала поиска решения и выбор численного метода реализовано с помощью диалогового окна.
Исходные данные могут задаваться как непосредственно, так и с помощью ссылок на ячейки электронной таблицы. Правая часть уравнения должна быть записана в соответствии с синтаксическими правилами MS Excel.
При запуске вычислительной процедуры осуществляется проверка корректности исходных данных (эта часть VBA кода не приведена в листинге ниже).
Полученное решение выводится в ячейки электронной таблицы (необходимо указать диапазон).
Листинг
часть программного кода кнопки "Решить"
'Инициализация
исходных данных
Points = (b - a) / step 'Вычисление количества точек
ReDim arr_x(Points) 'массив значений x
ReDim arr_y(Points) 'массив значений y в заданных точках
ReDim f(Points) 'массив значений производной в заданных точках
'ввод начальных
занчений х и у
arr_y(0) = y0
arr_x(0) = a
'вычисление
значения x с заданным шагом
For i = 1 To Points
arr_x(i) = arr_x(i - 1) + step
Next i
'Выбор метода
решения
If Opt_Euler.Value Then 'выбран метод Эйлера
For i = 1 To Points
arr_y(i) = euler(arr_x(i - 1), arr_y(i - 1), step, equat)
Next i
ElseIf Opt_Runge_kutta.Value Then 'выбран метод Рунге-Кутта
For i = 1 To Points
arr_y(i) = Runge_kutta(arr_x(i - 1), arr_y(i - 1), step, equat)
Next i
Else 'выбран метод прогноза и коррекции
'вычисление значений производной в начальной точке
f(0) = sheetFormula(arr_x(0), arr_y(0), equat)
'вычисление значения y в первых трех точках
'по методу Рунге-Кутта
For i = 1 To 3
arr_y(i) = Runge_kutta(arr_x(i - 1), arr_y(i - 1), step, equat)
f(i) = sheetFormula(CStr(arr_x(i)), CStr(arr_y(i)), equat)
Next i
'вычисление значения y в последующих точках
'по методу Милна
For i = 4 To Points
arr_y(i) = milna(arr_x(i), step, equat, i - 1, eps)
f(i) = sheetFormula(CStr(arr_x(i)), CStr(arr_y(i)), equat)
Next i
End If
'///////////////////////////////////////////////////////////////////
'Функция нахождения решения диффереренциального уравнения
' в заданной точке по методу Эйлера
'///////////////////////////////////////////////////////////////////
Function euler(ByVal x As Double, ByVal y As Double, ByVal step As
Double, ByVal eqaut As String) As Double
Dim y1 As Double, y0 As Double
f0 = sheetFormula(CStr(x), CStr(y), eqaut)
y1 = y + step * f0
euler = y + 1 / 2 * step * (f0 + sheetFormula(CStr(x + step), CStr(y1),
eqaut))
End Function
'///////////////////////////////////////////////////////////////////
'Функция нахождения решения диффереренциального уравнения
' в заданной точке по методу Рунге-Кутта
'///////////////////////////////////////////////////////////////////
Function Runge_kutta(ByVal x As Double, ByVal y As Double, ByVal step
As Double, ByVal eqaut As String) As Double
Dim y1 As Double, y0 As Double, k1 As Double, k2 As Double, k3 As
Double, k4 As Double
k1 = step * sheetFormula(CStr(x), CStr(y), eqaut)
k2 = step * sheetFormula(CStr(x + step / 2), CStr(y + k1 / 2), eqaut)
k3 = step * sheetFormula(CStr(x + step / 2), CStr(y + k2 / 2), eqaut)
k4 = step * sheetFormula(CStr(x + step), CStr(y + k3), eqaut)
Runge_kutta = y + (k1 + 2 * k2 + 2 * k3 + k4) / 6
End Function
'///////////////////////////////////////////////////////////////////
'Функция нахождения решения диффереренциального уравнения
' в заданной точке по методу Милна
'yp - значение Y, полученное по формуле прогноза
'yn1- значение Y, полученное с применением по формулы коррекции
'fp - производная по значению yp
'fn1 - производная по значению yn1
'eps - точность
'///////////////////////////////////////////////////////////////////
Function milna(ByVal x As Double, ByVal step As Double, ByVal eqaut As
String, ByVal n As Byte, ByVal eps As Double) As Double
Dim yp As Double
yp = arr_y(n - 3) + 3 / 4 * step * (2 * f(n) - f(n - 1) + 2 * f(n - 2))
fp = sheetFormula(CStr(x), CStr(yp), eqaut)
yn1 = arr_y(2) + step / 3 * (fp + 4 * f(n) + f(n - 1))
fn1 = sheetFormula(CStr(x), CStr(yn1), eqaut)
'вычисления выполняются до тех пор, пока не будет достигнута требуемая
точность
Do While Abs(fp - fn1) > eps
yp = yn1
fp = fn1
yn1 = arr_y(n - 1) + step / 3 * (fp + 4 * f(n) + f(n - 1))
fn1 = sheetFormula(CStr(x), CStr(yn1), eqaut)
Loop
milna = yn1
End Function
Пример
:
