Решение нелинейных уравнений в delphi

1. Делим интервал a..b пополам. Получаем c=(a+b)/2; Вычисляем f(c); Знак f(c) совпадает либо со знаком f(a) либо со знаком f(b) (либо f(c)=0 и с является корнем);
2. Если знак f(c) совпадает со знаком f(a), то корень находится на интервале c..b. Видим, что ситуация повторяется, поэтому делим пополам интервал c..b, то есть повторяем пункт a.
3. В противном случае, то есть если знак f(c) совпадает со знаком f(b), то корень находится на интервале a..c. Тогда делим пополам интервал a..b, то есть повторяем пункт a.
4. В результате интервал, в котором находится корень, уменьшается на каждом шаге в 2 раза, и это повторяется до тех пор, пока интервал не станет меньше заданной допустимой погрешности для значения корня.

• Запускаем среду Делфи7 (т е файл .\Program Files\Borland\Delphi7\Bin\delphi32.exe)
• Появляется несколько окон. Основные из них:

  Главное окно. Здесь имеется меню, панели с кнопками, панель компонентов (несколько закладок).
  Именно с панели компонентов мы перетаскиваем компоненты на форму. При запуске Делфи появляется пустая форма. Из пунктов меню часто используется File — для создания нового проекта (можно выбать тип проекта, в том числе — консольный), для открытия существующего проекта (при продолжении работы над проектом), для сохранения проекта.

  Это окно формы. Его размеры можно менять мышкой. Свойства формы, которые часто настраивают:
  — Caption (заголовок). По умолчанию — Form1. Можно (при помощи инспектора объектов — о нем дальше) заменить заголовок на «Это мой первый проект» или «Решение уравнений методом дихотомии» и т п.
  — Color т е цвет,
  — BorderStyle — Можно задать форму постоянных размеров (bsSingle) — Font.CharSet = RUSSIAN_CHARSET Отметим, что Font (шрифт), являющийся свойством формы, является объектом и мы здесь настроили свойство CharSet объекта Font. Часто определяют еще размер (Size) этого же объекта Font. Заметим, что настроенные свойства Font НАСЛЕДУЮТСЯ установленными на форму объектами.
  — Position иногда полезно, чтобы форма располагалась по центру экрана при запуске приложения. Для этого выберите Position = poScreenCenter из выпадающего списка возможных значений этого свойства.
  — TransparentColor, TransparentColorValue используйте их, если нужно часть формы сделать прозрачной, например, сделать круглую форму. Но при этом имейте в виду, что старые оп. системы (Windows 95,98 и т п) не поддерживают прозрачность.

  Это окно инспектора объектов. С его помощью настраиваются свойства формы и компонентов, установленных на форму. Вы можете выделить любой компонент на форме (с помощью мышиного щелчка), и свойства выделенного объекта будут видны и доступны для изменения в инспекторе объектов. В частности, на этой картинке инспектор объектов отображает свойства формы Form1. Видно, как программист установил свойство Position. Отметим также, что инспектор объектов имеет 2 закладки: Свойства (Properties) и События (Events). События — что-то новое для нас, поговорим о них немного.
  Примерами событий могут служить: щелчок мышкой, сбой при чтении файла или выполнении действий над данными, изменение позиции и размера окон и т п. События обнаруживаются операционной системой, которая, в свою очередь, информацию о событии посылает приложению. Компоненты Делфи способны реагировать на некоторые из событий. Список таких событий Вы увидите раскрыв закладку Events. Пока что мы будем пользоваться только событием onClick (мышиный щелчок по компоненту — кнопке). Говоря коротко, если сделать двойной щелчок в инспекторе — правее имени события, то в редакторе текста появится заготовка процедуры, в которой мы должны записать действия программы при наступлении этого события. Но об этом заботиться рано. Займемся пока текстом проекта.

  Делфи-проект состоит из модулей, и когда при запуске Делфи создается проект, то в окне редактора появляется заготовка текста модуля. Что интересного можно здесь увидеть?
  — Текст начинается с заголовка модуля unit Unit1; Видим что имя модуля (как и многие другие имена) создается автоматически (но мы можем переименовать).
  interface — объявляется начало интерфесной части модуля, (т е все, что здесь описано, можно сделать доступным в других модулях проекта.)
  Далее записано предложение uses: Здесь перечисляются модули, используемые в данном модуле. Не будем пока это обсуждать. Затем идут описания типов, в которых уже вписано описание КЛАССА TForm1 — это класс, описывающий создаваемую нами форму. Он будет наполняться при добавлении компонентов на форму. Думаю, что Вам полезно вспомнить, что такое класс.
  Далее следуют описания переменных, имеется всего 1 переменная: Form1: TForm; Это переменная — объект. Можно сказать, что типом для объекта является класс, так как описание деталек объекта содержит именно класс.

• Начнем добавление компонентов на форму. Посмотрите готовый проект. Как Вы видите, в левом верхнем углу формы находится текст:
  1. Если функция.
  Для этого перетащим компонент Label, установим его свойства: AutoSize = false, WordWrap = true, впишем текст в Caption и мышкой установим нужные размеры текста. Так как при настройке свойств формы было установлено: Font.Size = 10, то шрифт у Label1 тоже имеет Size = 10.
• Добавим рамку с заголовком «Задайте интервал a..b . » (компонент GroupBox1), зададим ей размеры и текст.
• Далее разместим внутри GroupBox1 компоненты Label2..5 для надписей a= b= f(a)= f(b)= а также Label6,7 — для показа значений f(a) f(b). В свойство Caption можно вписать 0.
• Нужно также поставить 2 компонента (класс TEdit) для ввода знаений a и b. Изменим имена этих редакторов на edA и edB (Это не обязательно, но удобно). Впишем свойство edA.Text = -0,2 и edB.Text = 2
• Теперь можно поставить кнопку (Button) с надписью «Построить график. «. Притащите ее с панели компонентов, впишите надпись в свойство Caption. Чтобы не путаться, присвоим кнопке имя (свойство name) btnCalc. Создадим также процедуру — обработчик щелчка по этой кнопке. Для этого: выделим кнопку, в инспекторе объектов перейдем на вкладку Events, найдем событие OnClick и сделаем двойной щелчок в правом окошке строки OnClick.
• Тогда в окне редактора появится шаблон-заготовка процедуры: Очень важно: если нужно удалить процедуру обработчика события, удалите из нее все, что Вы вписали, НО ОСТАВЬТЕ первоначальный шаблон. После этого просто сохраните проект и обработчики с пустым телом (т е пусто между begin и end) будут удалены автоматически. Если же Вы вручную удалите весь текст процедуры, то возникнет ошибка.
  Для начала впишем сюда действия по вычислению f(a) и f(b):
  Этот модуль может только вычислять функцию в двух заданных точках. В дальнейшем мы добавим все остальное. Но этот модуль работает и можно видеть результаты. Теперь нужно сохранить проект. Советую каждый проект сохранять в отдельной папке. Щелкните кнопку «Save All» или выберите в меню: File | Save All. Откроется окно, в котором выберите или создайте папку для проекта. Обратите внимание, что требуется 2 сохранения: для модуля и для проекта. При этом модуль и проект можно переименовать. После сохранения полезно заглянуть в эту папку и убедиться, что там записаны все файлы проекта (по умолчанию это Unit1.pas, Unit1.dfm, Project1.cfg, Project1.dof, Project1.dpr, Project1.res). Если эти файлы сохранены, то для продолжения работы с проектом необходимо и достаточно открыть Project1.dpr. Основное внимание мы уделяли файлу Unit1.pas — в нем содержится текст модуля Unit1. В тексте содержится описание класса TForm1, подпрограмма ff. Часть информации сохраняется в файле Unit1.dfm и именно поэтому нельзя удалять обработчики событий простым редактированием файла Unit1.pas — так как получается противоречие между Unit1.pas и Unit1.dfm. Остальные файлы проекта не нуждаются во внимании программиста, однако при переносе на другой компьютер могут «потеряться» файлы. В этом случае загляните в текст файла Project1.dpr и обратите внимание на предложение вроде:
  Unit1 in ‘Unit1.pas’ ;
  Здесь после in записывается путь к файлу Unit1.pas. Если все файлы проекта хранятся в одной папке, то путь к файлу пустой — указывается только имя файла. Если же Вы сохраняли проект не совсем правильно, то может записаться путь вроде: Unit1 in ‘D:\temp\Unit1.pas’ ;
  а на другой машине такой папки нет. Тогда подредактируйте Project1.dpr и укажите правильный путь. Лучше всего поместить все файлы проекта в одну папку и убрать пути.
  Проект можно компилировать (нажать Ctrl-F9). При этом могут появиться сообщения об ошибках.

  После устранения ошибок компиляция завершится нормально и в папке появится исполняемый файл Project1.exe. Этот файл можно запускать (открывать) на любом другом компьютере, если он работает под Windows (т е среда Делфи для этого не нужна). Исполняемые файлы не совсем простых программ могут не запускаться на другом компьютере, так как при компановке не подключены некоторые библиотеки. Эта проблема решается в Меню | Project | Options | Packages | Runtime Packages Также при компиляции не совсем простых программ могут «не найтись» некоторые библиотеки, не входящие в комплект Делфи. Такие библиотеки подключаются в Меню | Project | Options | Directories | Search Path — здесь добавляются в проект пути к дополнительным библиотекам.

  Можно совместить компиляцию с запуском на выполнение — нажмите клавишу F9.
  Итак, проект запущен и работает. Теперь можно подбором найти такие близкие значения a и b, чтобы f(a) и f(b) имели разные знаки + — и затем применить метод дихотомии на этом интервале a..b. Именно этим мы сейчас и займемся, а построение графика сделаем позже.

4.8 Решение систем нелинейных уравнений стр.1

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

4.8.1 Методы решения

В разд. 4.4 и 4.5 были описаны одномерные варианты методов простых итераций и Ньютона. В данном разделе рассмотрим многомерные варианты этих методов, когда вместо одного уравнения и одной переменной х имеется система п уравнений и вектор переменных

Метод простых итераций в многомерном случае применяется к системе уравнений, заданных в форме (4.2). Итерации организуются так же, как в одномерном случае:

Только, в отличие от одномерного варианта, X — это вектор переменных, §¥л — вектор функций.

В качестве критерия окончания используется оценка степени приближения к корню за итерацию: какая-либо норма ||Х^к+1) — Х )|тах. Если |^(Х(к+1>)|тах > |^(Х(к))|тах, значит метод начинает расходиться. Тогда параметр t уменьшается (например, в 2 раза) и происходит возврат на шаг 3.

Шаг 6. Проверяются критерии окончания. Если они выполнены, то расчет завершается. В противном случае переходом на шаг 1 начинается следующая итерация.

Как видно из приведенного алгоритма, если невязка все время уменьшается, то параметр Ь всегда равен 1 и алгоритм работает как обычный метод Ньютона. Но если невязка возросла, то повторяются шаги 3-5 со все уменьшающимся параметром пока невязка на данной итерации не станет меньше или хотя бы равной невязке предыдущей итерации. Таким образом, данный метод не может расходиться, так как невязка просто не может увеличиваться от итерации к итерации.

Повторение шагов 3-5 при возникновении опасности расходимости — не слишком трудоемкий процесс. Он не связан с вычисление поправки Ньютона А и, значит, не требует расчетов производных и решения системы (4.14). Однако так как теоретически цикл по шагам 3-5 может оказаться бесконечным, надо ограничить минимальное значение параметра Если t станет меньше некоторого предельного значения, следует прервать расчет, известив пользователя об отсутствии сходимости.

Об автоматическом дифференцировании, методе Ньютона и решении СЛАУ на 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).