Численные методы решения дифференциальных уравнений в delphi

Об автоматическом дифференцировании, методе Ньютона и решении СЛАУ на Delphi. Часть 1

Об автоматическом дифференцировании (АД) на Хабре уже писалось здесь и здесь. В данной статье предлагается реализация АД для Delphi (протестировано в Embarcadero XE2, XE6) вместе с удобными классами методов Ньютона для решения нелинейных уравнений f(x) = 0 и систем F(X) = 0. Любые ссылки на готовые аналогичные библиотеки приветствуются, сам же я подобного не нашел, не считая отличного решателя СЛАУ с разреженной матрицей (см. под катом).

Давайте в самом начале отметим для себя, что выбор Delphi обусловлен legacy-кодом, тем не менее на C++ задачу можно решать следующим образом. Во-первых, для методов Ньютоновского (базовый метод Ньютона, метод Бройдена) типа имеются Numerical Recipes in C++. Ранее «Рецепты» были только на чистом C и приходилось делать свои классовые обертки. Во-вторых, можно взять одну из АД-библиотек из списка Autodiff.org. По моему опыту использования CPPAD быстрее FADBAD и Trilinos::Sacado примерно на 30%, самая же быстрая, судя по описанию, библиотека это новая ADEPT. В-третьих, для матрично-векторных операций можно взять проверенный временем uBlas , либо новые и быстрые конкуренты Armadillo и blaze-lib — это если для решения СЛАУ использовать отдельные библиотеки (например, SuiteSparce или Pardiso для прямых и ITL для итерационных методов). Привлекательным является также использование интегрированных библиотек линейной алгебры и решателей СЛАУ Eigen, MTL, PETSc (имеются также Ньютоновские решатели на C). Вся триада из заголовка полностью реализована, насколько мне известно, только в Trilinos.

О применении численного дифференцирования

Производные можно вычислять аналитически и численно. К аналитическим методам относятся ручное дифференцирование, символьное (Maple, Wolfram и т.п.) и непосредственно автоматическое дифференцирование, выраженное в средствах выбранного языка программирования.

Современный тренд на использование АД оправдан одной простой причиной — с помощью этой техники устраняется избыточность кода и его дублирование. Другой аргумент состоит в том, что, например, при решении нелинейных дифференциальных уравнений (систем) сеточными методами способ вычисления F(X) сам по себе является нетривиальной задачей. В реальных задачах невязка F(X) представлена суперпозицией вызовов функций с разных слоев программы и ручное дифференцирование теряет свою наглядность. Следует также отметить, что при моделировании нестационарных процессов F(X) меняется на каждом шаге по времени, также может меняться и сам вектор неизвестных X. Использование АД позволяет нам сконцентрироваться непосредственно на формировании F(X), однако не снимает вопрос о верификации получаемой матрицы Якобиана dF(X)/dX. Дело в том, что при вычислении невязок мы можем забыть изменить тип промежуточной переменной со стандартного double на тип АД (а многие библиотеки имеют неявное приведение типа АД к double), тем самым некоторые производные будут вычислены некорректно. Проблема в этом случае состоит в том, что даже при наличии ошибок в формулах для производных метод Ньютона может сходиться, хоть и за возросшее число итераций. Это может быть незаметно при одних начальных данных и приводить к расходимости процесса при других.

Таким образом, какой бы аналитической способ дифференцирования df/dx не был выбран, его крайне желательно дополнить сравнением с численным дифференцированием (f(x + h) — f(x)) / h, иначе всегда будут оставаться сомнения в правильности кода. Естественно, численные производные никогда не совпадут с точностью с правильными аналитическими, тем не менее можно порекомендовать следующий прием юнит-тестирования. Можно выгрузить в текстовые файлы вектора X, F(X) и матрицу dF(X)/dX и выложить на SVN. Затем получить dF(X)/dX численно и сохранить файл поверх существующего, после чего визуально сравнивать файлы между собой. Здесь Вы всегда увидите насколько поменялись значения и сможете локализовать координаты элементов матрицы с большими отклонениями (не в долях) — в этом месте находится ошибка аналитического дифференцирования.

Реализация АД

В Embarcadero Delphi до версии XE5 отсутствует возможность перегрузки арифметических операций для классов, но есть такая возможность для структур record (ссылка):

Обычно в АД на C++ размерность вектора производных является переменной величиной и инициализируется в конструкторе. В Delphi нельзя (есть попытки обойти) перегрузить оператор присваивания для структур и в связи с побитовым копированием вектор производных приходится делать фиксированной длины. Соответствующая константа TAutoDiff.n_all должна изначально задаваться вручную.

Пример 1

На данный момент в библиотеке перегружены почти все бинарные и унарные операторы, за исключением тригонометрических функций и булевых операций сдвига. Недостающие функции легко доработать самостоятельно.

Реализация АД для разреженных матриц

С одной стороны фиксированное значение n_all это существенное ограничение, ведь размерность вектора может поступать извне. С другой стороны для некоторых задач его можно ослабить. Дело в том, что если говорить о численных методах решения уравнений механики сплошных сред, то возникающие в них матрицы имеют разреженную структуру — классический пример это «схема крест» для оператора Лапласа, когда в уравнении для узла с координатами (i, j) (ограничимся 2D) задействованы только 5 узлов: (i, j), (i-1, j), (i+1, j), (i, j-1), (i, j+1). Обобщая идею мы должны заложить следующее для данной конкретной задачи:

Пусть общее число узлов в решаемой задаче N. Тогда в матрице Якобиана N_all = N * n_junc_vars столбцов, из них ненулевых только n_all. Если завести теперь внутри структуры TAutoDiff целочисленный вектор n_juncs, каждый элемент которого определяет глобальный индекс узла от 0 до N-1, то функцию dx с локальным индексом из диапазона [0, n_all-1] можно дополнить функцией dx_global с уже глобальным индексом из диапазона [0, N_all-1]. Пример 2 иллюстрирует детали использования такого интерфейса, плюсы которого будут наиболее полно видны при реализации метода Ньютона ниже.

Пример 2

В следующей части планируется рассмотрение класса методов Ньютоновского типа, а также вопроса выбора разреженного решателя СЛАУ. Пока же предлагаю читателям:

  • попробовать написать АД на C++11 с использованием семантики перемещений: 1) это должно работать очень быстро; 2) можно будет обойтись перегрузкой операторов без expression templates; 3) это будет впервые.
  • идею для курсовой по реализации АД на Roslyn CTP: можно работать сразу с синтаксическим деревом, которое содержит всю информацию об арифметических выражениях в F(X).

Метод Рунге-Кутта решения дифференциальных уравнений и их систем

Delphi site: daily Delphi-news, documentation, articles, review, interview, computer humor.

Метод Рунге-Кутта решения дифференциальных уравнений и их систем

Метод позволяет решать системы обыкновенных дифференциальных уравнений (ОДУ) первого порядка следующего вида:
,
,
и т.д.,

которые имеют решение:

,
,
и т.д.,

где t – независимая переменная (например, время); X, Y и т.д. – искомые функции (зависимые от t переменные). Функции f, g и т.д. – заданы. Также предполагаются заданными и начальные условия, т.е. значения искомых функций в начальный момент.

Одно дифференциальное уравнение – частный случай системы с одним элементом. Поэтому, далее речь пойдет для определенности о системе уравнений.

Метод может быть полезен и для решения дифференциальных уравнений высшего (второго и т.д.) порядка, т.к. они могут быть представлены системой дифференциальных уравнений первого порядка.

Метод Рунге-Кутта заключается в рекурентном применении следующих формул:


.
где
,
,
,
,
,
,
,

Реализация Метода Рунге-Кутта на Delphi может выглядеть так (привожу полностью модуль):

type
TVarsArray = array of Extended; // вектор переменных включая независимую
TInitArray = array of Extended; // вектор начальных значений
TFunArray = array of function(VarsArray: TVarsArray ):Extended;
// вектор функций
TResArray = array of array of Extended; // матрица результатов
TCoefsArray = array of Extended; // вектор коэффициетов метода

function Runge_Kutt( // метод Рунге-Кутта
FunArray: TFunArray; // массив функций
First: Extended; // начальная точка по независимой координате
Last: Extended; // конечная точка по независимой координате
Steps: Integer; // число шагов по независимой координате
InitArray: TInitArray; // вектор начальных значений
var Res: TResArray // матрица результатов включая независ. переменную
):Word;
// возвращаемое значение — код ошибки

implementation
Function Runge_Kutt( // метод Рунге-Кутта
FunArray: TFunArray; // массив функций
First: Extended; // начальная точка по независимой координате
Last: Extended; // конечная точка по независимой координате
Steps: Integer; // число шагов по независимой координате
InitArray: TInitArray; // вектор начальных значений
var Res: TResArray // матрица результатов включая независ. переменную
):Word; // возвращаемое значение — код ошибки
var
Num: Word; // число уравнений
NumInit: Word; // число начальных условий
Delt: Extended; // шаг разбиения
Vars: TVarsArray; // вектор переменных включая независимую
Vars2,Vars3,Vars4: TVarsArray; // значения перем. для 2-4 коэф.
Coefs1: TCoefsArray; // вектор 1-ыx коэффициентов в методе
Coefs2: TCoefsArray; // вектор 2 коэффициентов в методе
Coefs3: TCoefsArray; // вектор 3 коэффициентов в методе
Coefs4: TCoefsArray; // вектор 4 коэффициентов в методе
I: Integer; // счетчик цикла по иттерациям
J: Word; // индекс коэф.-тов метода
K: Integer; // счетчик прочих циклов
begin
Num:=Length(FunArray); // узнаем число уравнений
NumInit:=Length(InitArray); // узнаем число начальных условий
If NumInit<>Num then
begin
Result:=100; // код ошибки 100: число уравнений не равно числу нач. усл.
Exit;
end;
Delt:=(Last-First)/Steps; // находим величину шага разбиений
SetLength(Res,Num+1,Steps+1); // задаем размер матрицы ответов с незав. перем.
SetLength(Vars,Num+1); // число переменных включая независимую
SetLength(Vars2,Num+1); // число переменных для 2-го коэф. включая независимую
SetLength(Vars3,Num+1); // число переменных для 3-го коэф. включая независимую
SetLength(Vars4,Num+1); // число переменных для 4-го коэф. включая независимую
SetLength(Coefs1,Num); // число 1-ыx коэф. метода по числу уравнений
SetLength(Coefs2,Num); // число 2-ыx коэф. метода по числу уравнений
SetLength(Coefs3,Num); // число 3-иx коэф. метода по числу уравнений
SetLength(Coefs4,Num); // число 4-ыx коэф. метода по числу уравнений
// Начальные значения переменных:
Vars[0]:=First;
For K:=0 to NumInit-1 do Vars[K+1]:=InitArray[K];
For J:=0 to Num do Res[J,0]:=Vars[J]; // первая точка результата
For I:=0 to Steps-1 do // начало цикла иттераций
begin
For J:=0 to Num-1 do Coefs1[J]:=FunArray[J](Vars)*delt; // 1-й коэфф.
// Находим значения переменных для второго коэф.
Vars2[0]:=Vars[0]+delt/2;
For K:=1 to Num do Vars2[K]:=Vars[K]+Coefs1[K-1]/2;
For J:=0 to Num-1 do Coefs2[J]:=FunArray[J](Vars2)*delt; // 2-й коэф.
// Находим значения переменных для третьго коэф.
Vars3[0]:=Vars[0]+delt/2;
For K:=1 to Num do Vars3[K]:=Vars[K]+Coefs2[K-1]/2;
For J:=0 to Num-1 do Coefs3[J]:=FunArray[J](Vars3)*delt; // 3 коэфф.
// Находим значения переменных для 4 коэф.
Vars4[0]:=Vars[0]+delt;
For K:=1 to Num do Vars4[K]:=Vars[K]+Coefs3[K-1];
For J:=0 to Num-1 do Coefs4[J]:=FunArray[J](Vars4)*delt; // 4 коэфф.
// Находим новые значения переменных включая независимую
Vars[0]:=Vars[0]+delt;
For K:=1 to Num do
Vars[K]:=Vars[K]+(1/6)*(Coefs1[K-1]+2*(Coefs2[K-1]+Coefs3[K-1])+Coefs4[K-1]);
// Результат иттерации:
For J:=0 to Num do Res[J,I+1]:=Vars[J];
end; // конец итераций
Result:=0; // код ошибки 0 — нет ошибок
end;

Модуль полностью работоспособен. Возвращаемое функцией Runge_Kutt значение – код ошибки. Вы можете дополнить список ошибок по своему усмотрению. Рассчитанные функции системы помещаются в массив Res. Чтобы не загромождать код, в модуле опущены проверки (типа блоков try). Рекомендую их добавить по своему усмотрению.

Ниже приводится описание функции Runge_Kutt и типов, использующихся в модуле.

Function Runge_Kutt (FunArray: TFunArray; First: Extended; Last: Extended; Steps: Integer; InitArray: TInitArray; var Res: TResArray):Word;

Здесь:
FunArray — вектор функций (правых частей уравнений системы);
First, Last — начальная и конечная точки расчетного интервала;
Steps — число шагов по расчетному интервалу;
InitArray — вектор начальных значений
Res — матрица результатов включая независимую переменную.

В модуле описаны типы:

type
TVarsArray = array of Extended; // вектор переменных включая независимую
TInitArray = array of Extended; // вектор начальных значений
TFunArray = array of function(VarsArray: TVarsArray ):Extended; // вектор функций
TResArray = array of array of Extended; // матрица результатов
TCoefsArray = array of Extended; // вектор коэффициетов метода

Функция возвращает коды ошибок:
0 – нет ошибок;
100 — число уравнений не равно числу начальных условий.

Решение содержится в переменной-матрице Res. Первый индекс матрицы относится к переменной (0 – независимая переменная, 1 – первая зависимая и т.д.), второй – к номеру расчетной точки (0 – начальная точка).

Рассмотрим один пример использования модуля. Создадим новое приложение и подключим к нему модуль. На форме приложения разместим кнопку Button1 и область текста Memo1. Поместим в приложение две функции и обработчик нажатия кнопки:

//Задаем функции (правые части уравнений)
function f0(VarArray:TVarsArray):extended;
begin
Result:=4*VarArray[0]*VarArray[0]*VarArray[0];
end;

function f1(VarArray:TVarsArray):extended;
begin
Result:=1;
end;

procedure TForm1.Button1Click(Sender: TObject);
var
I: Integer;
FunArray: TFunArray; // массив функций
First: Extended; // начальная точка по независимой координате
Last: Extended; // конечная точка по независимой координате
Steps: Integer; // число шагов по независимой координате
InitArray: TInitArray; // вектор начальных значений
Res: TResArray; // матрица результатов включая независ. переменную
begin // Создаем вектор функций:
SetLength(FunArray,2);
FunArray[0]:=f0;
FunArray[1]:=f1;
// Задаем интервал и число шагов:
First:=0;
Last:=10;
Steps:=10;
// Задаем начальные условия:
SetLength(InitArray,2);
InitArray[0]:=0;
InitArray[1]:=0;
// Вызов метода и получение результатов:
Memo1.Lines.Clear;
I:=Runge_Kutt(FunArray, First, Last, Steps, InitArray, Res);
ShowMessage(‘Код ошибки = ‘+IntToStr(I));
For I:=0 to Steps do
Memo1.Lines.Add(floattostr(Res[0,I])+’ ‘+floattostr(Res[1,I])+’ ‘+floattostr(Res[2,I]));
end;

Модуль с примером и справкой можно скачать: русский вариант (бесплатно) или английский вариант (условно-бесплатно).

Copyright© 2006 Андрей Садовой Специально для Delphi Plus

Численное решение задачи Коши (ЗК) для системы обыкновенных дифференциальных уравнений (Delphi)

Работа сделанна в 2004 году

Численное решение задачи Коши (ЗК) для системы обыкновенных дифференциальных уравнений (Delphi) — Курсовая Работа, раздел Программирование, — 2004 год — Санкт-Петербургский Государственный Университет Факультет Прикладной Математи.

Санкт-Петербургский Государственный Университет Факультет Прикладной Математики — Процессов Управления Курсовая работа по Методам Вычислений на тему Численное решение задачи Коши ЗК для системы обыкновенных дифференциальных уравнений Выполнил Студент 2 курса 20 группы Орлов Евгений Вадимович Проверил Санкт-Петербург 2004г. ОГЛАВЛЕНИЕ Явный метод типа Рунге-Кутта. 4 Многошаговые методы. 5 Правило Рунге для практической оценки ГП. 5 Алгоритм автоматического выбора шага. 5 Текст программы. 7 Выходные данные 10 Явный метод типа Рунге-Кутта. 1 где — искомая, а — заданная векрот-функция, -заданный вектор. заданна и непрерывна в области , Найти решение системы, определнной на конечном отрезке и удовлетворяющее начальному условию 2. Явные одношаговые методы типа Рунге-Кутта численные методы, в которых приближнное решение задачи 1,2 вычисляется по формулам Для метода пятого порядка точности S6, Остальные параметры этого метода совпадают с одноименными параметрами метода 4 порядка точности. узлы сетки шаг сетки. i2 i3 i4 i5 i6 Величина называется погрешностью метода в узле или глобальной погрешностью метода.

Если, то сетка называется равномерной.

Многошаговые методы

Многошаговые методы. В многошаговых методах обычно используют равномерную сетку, так как расчтные формулы в этом случае значительно проще. Пусть. Экстаполяционный метод Адамса явный. K шаговый ЭМА имеет порядок точности k. Формула линейного k1 шагового экстаполяционного метода Адамса имеет вид Таким образом, для четырехшагового ЭМА формула примет вид , Порядок точности четырехшагового ЭМА 4. Величина называется погрешностью метода в узле или глобальной погрешностью метода.

Правило Рунге для практической оценки ГП

Вычислим одним и тем же численным методом два приближнных решения зада. При реализации численных методов решения задачи Коши погрешность прибл. Если при изменении алгоритма автоматического выбора шага в многошагово. Исходная сетка сгущается до тех пор, пока не станет где — заданная гра. Согласно Таким образом.

Алгоритм автоматического выбора шага

Алгоритм автоматического выбора шага. б число узлов должно быть наименьшим. В противном случае делится пополам, и с этим новым значением шага все . точка с первым значением шага, при котором выполняется неравенство, бе. Пусть узлы построены. При построении узла учтем требование Б если, то .

Текст программы

4h1-1.2-3dyy0-8.5yy0-1.2 for i 3 to p2 do rezult31,irezult21,i-1h11223. Текст программы. APPTYPE CONSOLE F чтобы можно было передавать функции. 5yyi-1h13.4i-1h1-1.2-163dyyi-2h1-8.5yyi- 2h13.4i-2h1-1.253dyyi-3h1-8.5.

Выходные данные

Hr160 Runge t 0.1, y1t 0.9893166478, y2t 1.2099825961, t 0.2, y1t 0.95. Выходные данные.


источники:

http://www.delphiplus.org/articles/algorithm/rk_method/index.html

http://allrefs.net/c20/1bax2/