Введение. 
Компьютерные технологии (программирование на C#)

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

В настоящее время на смену прежней библиотеке функций и структур Win32, созданной в начале 90-ых годов и составляющей основу OS Windows, пришла новая библиотека .NET (читается «дот-нет»; можно перевести как «сетевое расширение»).

Библиотека .NET создана в начале 2000;ых годов и является библиотекой классов. В настоящее время библиотека .NET составляет основу операционной системы нового поколения. Настоящий курс знакомит с библиотекой .NET и языком программирования C# (читается «си-шарп»; можно перевести как «на пол тона выше си», или «си-диез»), лежащим в ее основе.

Неформальное введение в C#. Интерфейс IIntegrator.

При численной симуляции физических процессов часто встречается ситуация, в которой математическая модель предполагает решение задачи с начальными условиями, или задачи Коши. Поэтому, в качестве первого примера, рассмотрим, как может выглядеть класс интегрирования системы дифференциальных уравнений на языке C#.

Формулировка задачи Коши выглядит довольно просто.

Дана система N обыкновенных дифференциальных уравнений 1-ого порядка вида.

dyi/dt = fi (t, y1, y2,…, yN).

с N неизвестными функциями yi (t).

Даны начальные условия — значения неизвестных функций yi0 в некоторой точке t = t0.

Необходимо определить значения функций yi в любой точке t.

Правые части дифференциальных уравнений определены задачей. В каждой точке с координатами t, y1, y2,…, yN функции fi образуют вектор относительно линейных преобразований переменных dyi. Эти вектора fi (t, y1, y2,…, yN), заданные в каждой точке, образуют, как говорят, векторное поле.

При составлении алгоритма необходимо предусмотреть, чтобы класс интегратора имел доступ к следующим данным Число уравнений N;

Метод, позволяющий вычислять векторное поле fi (t, y1, y2,…, yN) в каждой точке;

Начальные условия;

Погрешность, с которой проводится интегрирование;

Метод, реализующий алгоритм интегрирования до заданной точки;

Наконец, метод, который предоставит пользователю возможность отладки и управления кодом в ходе процесса интегрирования. Это так называемая «функция обратного вызова» (callback function), или «обработчик события».

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

Язык C# позволяет уже на этом, столь общем этапе программирования задачи, написать классы, которые позволят ограничить рамки задачи и кода, требуемого для ее решения. На основании сформулированных требований можно составить так называемый «интерфейс» — совершенно абстрактный класс, состоящий лишь из заголовков методов, которые должны быть реализованы в любом классе, решающем поставленную задачу.

Назовем наш интерфейс IIntegrator. В языке C# все имена интерфейсов принято начинать с буквы I. Синтаксис описания интерфейса выглядит следующим образом.

public interface IIntegrator.

{.

// Здесь описываются члены интерфейса.

}.

Служебные слова в дальнейшем выделяются полужирным шрифтом. Описание интерфейса IIntegrator вместе с некоторыми другими вспомогательными классами будет нами в дальнейшем помещено в отдельную библиотеку классов. Служебный термин public является одним из так называемых модификаторов доступа. Он означает, что описание интерфейса будет доступно любому приложению, ссылающемуся на ту библиотеку, в которой это описание размещается.

Описание так называемых членов интерфейса должно размещаться внутри фигурных скобок, которые ограничивают описание интерфейса от остальных классов, описываемых в том же модуле (библиотеке). Любой комментарий к коду, который не должен влиять на его трансляцию, может быть помещен после двойного знака // (double slash) и до конца строки.

В интерфейсе IIntegrator должно быть реализовано свойство — число уравнений N. Этот факт мы опишем в качестве одного из членов интерфейса — так называемого свойства (property).

ushort N { get;}.

Так в C# выглядит описание абстрактного (не реализованного) свойства.

Все члены интерфейсов имеют доступ public, и этот модификатор нигде в описании членов специально не указывается.

Свойства в C# - это пара методов доступа set и get. Методы доступа устанавливают (set) и возвращают (get) объект того типа, который указан в описании свойства. Метод установки set имеет неявный параметр, обозначаемый служебным словом value. Параметр value является объектом типа свойства. Он содержит в себе значение, которое передается свойству оператором присваивания, если имя свойства находится в правой части оператора. В языке C# оператор присваивания имеет вид простого знака равенства =. Поэтому запись N = 5, могла бы означать в коде C#, что некоторому свойству с именем присваивается значение 5.

Сразу отметим, что в нашем случае свойства N такое присваивание оказалось бы невозможным. В описании свойства N отсутствует метод set. А это означает, что интерфейсом IIntegrator предполагается реализация свойства N только по чтению. Метод get свойства возвращает объект типа свойства, то есть «работает» в операторе присваивания, где свойство находится в правой части. Например, curN = N, где curN предположительно есть некоторое текущее (current) значение свойства N.

Тип ushort описанного свойства N означает, что свойство N должно возвращать любое целое неотрицательное число, меньшее 65 536.

Приведенное описание свойства числа уравнений N в форме «свойства только для чтения» предполагает, что объекты классов, основанных на интерфейсе IIntegrator (говорят, «классы, наследующие интерфейс»), уже при своем конструировании будут иметь заданное число уравнений, для решения которых они создаются. Поэтому необходимости в изменении значения числа уравнений на протяжении жизненного цикла таких объектов не будет возникать.

То же замечание будет относиться и к следующему свойству интерфейса IIntegrator.

TEquations Equations { get;}.

Предполагается, что свойство Equations будет возвращать ссылку на метод, который готовит локальный вектор уравнений. Тип этого метода имеет имя TEquations и должен быть описан где-либо вне скобок, содержащих члены интерфейса IIntegrator.

Описание типа TEquations имеет вид.

public delegate void TEquations.

(.

double t, //независимая переменная.

double [ ] y, //массив текущих значений неизвестных функций.

double [ ] localVector //массив правых частей уравнений.

);

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

Служебное слово delegate означает, что тип TEquations принадлежит к специальным классам (делегатам), полями которых являются ссылки на методы заданной сигнатуры. Термин сигнатура означает совокупность параметров метода, их типы и тип возвращаемого методом результата. При описании класса делегата сигнатура указывается в самом описании. Так, в данном случае, полями делегата типа могут быть ссылки на любой метод, который ничего не возвращает. На это указывает термин void. Метод, на который ссылается делегат, должен иметь три параметра. Первый параметр t должен быть типа double, а два других y и localVector должны быть одномерными массивами объектов типа double произвольной длины. На это указывают скобки [ ]. Объекты типа double хранят положительные и отрицательные числа с плавающей запятой, модуль которых лежит в интервале [5.0Ч10−324; 1.7Ч10 308].

Из комментария ясен смысл этих параметров. Отметим, что параметр типа double передается методу по значению. Другими словами передаются поле структуры, описанной в библиотеке .NET под именем Double, содержащее значение параметра t и длиной 8 байт. В то же время параметры типа массива double [ ] передаются в метод по ссылке, то есть передается адрес области памяти, в которой расположены элементы массива. Длина адреса в системе с 32-битовыми адресами 4 байта.

Итак, свойства N и Equations описывают систему дифференциальных уравнений. Предполагается, что эти свойства будут возвращать определенные значения после создания объекта реализующего (наследующего) интерфейс IIntegrator. При этом свойство N будет возвращать именно значение, содержимое 2-байтового поля структуры UInt16 библиотеки .NET (ushort — псевдоним этого типа), а свойство Equations — 4-ех байтовую ссылку на экземпляр класса делегата типа TEquations.

Интерфейс IIntegrator может предложить для реализации свойства, которые используются для задания начальных условий интегрирования и чтения неизвестных функций при текущем значении аргумента после, или в процессе интегрирования уравнений. Все эти свойства можно объединить в одно свойство, снабженное индексом. Такое свойство называется индексатором. По правилам C# описание индексатора как члена интерфейса имеет вид.

double this [ushort index] { get; set;}.

Служебное слово this необходимо при описании индексатора. Оно означает, что имя свойства будет совпадать с именем объекта, свойством которого является. Конкретное значение параметра передается в качестве дополнительного параметра методу при установке значения интегратора, и методу при получении этого значения. Пусть, например, описан объект типа IIntegrator (или класса-наследника от IIntegrator) с именем integrator. Тогда оператор присваивания.

integrator[2] = -0.5;

вызовет метод set индексатора this, передав ему в качестве параметра index значение 2, а значение -0.5 в качестве параметра value. Так можно установить начальное значение независимой переменной, сопоставив ей индекс 0, и начальные значения неизвестных функций, сопоставив им индексы от 1 до N.

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

y3 = integrator[3];

В общем случае тип параметра index индексатора (в нашем случае ushort), как и тип возвращаемого объекта (в нашем примере double) могут быть любыми. Параметр index является не индексом элемента массива, а параметром методов set и get индексатора.

Интегрирование уравнений всегда производится с ограниченной погрешностью. Поэтому интерфейс IIntegrator предлагает для реализации член типа.

double Tolerance { get; set;}.

Класс-наследник должен реализовать оба метода доступа set и get свойства Tolerance (погрешность).

Следующее поле.

CallBackEventHandler CallBack { get; set;}.

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

Дело в том, что метод, реализующий алгоритм интегрирования дифференциальных уравнений, может требовать выполнения довольно большого цикла вычислений. Представьте себе, что хотя бы одна из неизвестных функций yi (t) имеет период существенного изменения значительно меньший интервала между начальной и конечной точкой интегрирования. Тогда, для поддержки необходимой погрешности, цикл вычислений будет состоять из большого количества шагов. Возможно, пользователю будет необходимо выводить текущую информацию где-либо внутри цикла интегрирования, или даже на каждом его шаге. В этом случае класс, реализующий метод интегрирования, может предусмотреть вызов пользователем так называемой функции обратного вызова (callback function). Поставив вызов этой функции в некоторые точки цикла вычислений (например, в конце каждого шага цикла), можно обеспечить пользователя доступом к этому событию (например, окончания шага цикла). Свойство CallBack устанавливает и возвращает ссылку на функцию обратного вызова. Эта функция является объектом класса-делегата, названного здесь CallBackEventHandler — обработчик события обратного вызова.

Описание класса CallBackEventHandler следует поместить вне скобок с членами интерфейса IIntegrator.

public delegate void.

CallBackEventHandler (Object Sender, CallBackEventArgs e);

Как и в случае класса-делегата TEquations в этом описании указана сигнатура методов, которые могут быть делегированы типом CallBackEventHandler. Согласно описанию у метода (или методов), который может делегировать пользователь в качестве обработчика рассматриваемого события, должно быть два параметра типа Object и CallBackEventArgs и этот метод не должен возвращать какой-либо объект (void). Такой набор параметров является стандартным в модели обработки событий платформы .NET, хотя и не обязательным.

Тип Object является библиотечным классом — предком всех других классов библиотеки .NET. Поэтому на место параметра Sender (посыльный) может подставляться объект любого типа. Посыльным является объект, порождающий событие, с которым связан вызов метода-делегата, или обработчика события. Таким объектом может быть integrator, описанный выше. Этот объект вызывает обработчик внутри своего цикла и передает ему в качестве параметра Sender ссылку на себя (this). Ниже мы увидим этот фрагмент кода.

Второй параметр обработчика, названный e (от слова event — событие), имеет тип CallBackEventArgs (аргументы события обратного вызова).

Описание класса CallBackEventArgs, которое мы разместили в той же библиотеке, что и весь предыдущий код, имеет вид.

public class CallBackEventArgs: EventArgs.

{.

public bool Stop;

}.

Это пример описания класса, который имеет предка. Имя класса-предка указывается после имени класса-наследника через двоеточие.

Здесь предком является класс EventArgs, взятый нами из библиотеки .NET. Объекты класса EventArgs и его наследников являются стандартными параметрами обработчиков событий библиотеки .NET. Эти классы содержат данные о событии. Класс-предок пуст в том смысле, что его объекты не имеют полей — абстрактное событие не имеет данных. Но классы-наследники могут содержать поля, переносящие данные о событии.

В нашем случае наследник CallBackEventArgs имеет одно поле Stop типа bool. Поле Stop может принимать два значения true и false. По умолчанию после создания объекта класса CallBackEventArgs значение поля Stop равно false. Но пользователь, написав собственный код обработчика, может изменить это значение на true. В реализации интерфейса IIntegrator значение true может служить сигналом к прерыванию цикла интегрирования. В этом случае пользователь сможет управлять результатом события, меняя поле stop параметра e метода, делегированного в качестве обработчика событий.

На этом моменте мы еще раз остановимся при обсуждении кода реализации интерфейса. Здесь отметим, что в сигнатуре метода, который может быть делегирован в качестве обработчика событий делегату CallBackEventHandler, типом второго параметра может быть не только класс CallBackEventArgs, но и его предок — класс EventArgs. Это допускается синтаксисом и называется свойством контравариантности делегата. Правда в этом случае пользователь лишит себя возможности использовать поле stop, которое отсутствует у класса EventArgs.

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

bool IntegrateTo (double tEnd);

Метод IntegrateTo проводит интегрирование до точки tEnd и возвращает true или false в зависимости от того, остановлено ли его выполнение обработчиком (false), или метод завершен естественным путем (true) и найдены значения неизвестных функций в точке tEnd.

Посмотрите прокомментированный код описанных классов на файле IIntegrator.txt.

В описании кода присутствует строка.

using System;

Это ссылка на то, что в файле будут использованы классы из пространства имен System. Это пространство имен одной из стандартных библиотек .NET. В частности, в пространстве имен System описан использованный нами класс EventArgs.

Пространства имен используются при написании кода для того, чтобы логически разделить участки кода. Например, отделить имена классов одной библиотеки от имен классов другой. В частности, весь написанный нами код помещен в пространство имен namespace Integrators. Это имя, которое мы дали создаваемой нами библиотеки интеграторов. Для внешнего пользователя имя интерфейса IIntegrator выглядит как Integrators.IIntegrator. Оператор точка соединяет в полное имя идентификатора класса IIntegrator и имя пространства имен Integrators. Если пользователь в каком-либо файле кода будет ссылаться на классы нашей библиотеки, то в начале файла он сможет набрать строку.

using Integrators;

А затем использовать краткие имена наших классов IIntegrator, CallBackEventHandler и т. д.

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

Например,.

/// Класс аргументов обработчика события, наступающего внутри.

/// цикла интегрирования (например, после завершения оптимального шага).

Этот комментарий используется средой программирования для поддержки справочной системы библиотеки создаваемых классов. Другими словами, давая такой комментарий, мы создаем основу для справочной системы своего кода.

Идентификаторы, окруженные угловыми скобками, называются «тэгами». Эти тэги являются составной частью универсального языка общения, названного XML (eXtensible Markup Language — «расширяемый язык разметок»). Язык XML широко используются средой для обмена данными между ее отдельными компонентами. В отличие от известного языка HTML, построенного на той же идее разметок-тэгов, язык XML использует любые тэги, заданные его схемой. Кроме того, тэги языка XML используется не для указания формата изображаемых объектов, а для определения содержательной части — данных. Так в нашем примере «открывающий» и «закрывающий» тэги говорят о том, что между ними находится краткая справка о классе, описанном ниже. Везде в коде при ссылке на класс, обладающий таким комментарием, среда будет выводить справку на экран.

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

/// Хранит флаг прерывания интегрирования.

/// Служит для прерывания интегрирования, если установлено в true.

Еще один пример

/// Интегрирование от текущей точки до tEnd.

/// Значение независимой переменной,.

/// в котором интегрирование должно быть прервано.

/// Обычное завершение — возвращает true,.

/// завершение, прерванное обработчиком внутреннего события — false.

показывает, как описание метода может сопровождаться справкой о его параметрах (тэг.

Реализация интегратора абстрактным классом Рассмотрим первый этап реализации интерфейса IIntegrator, спроектировав пока что абстрактный класс TIntegrator, наследующий этот интерфейс.

Описание класса синтаксически похоже на описание интерфейса. Класс, как и интерфейс, имеет заголовок с именем и модификатором доступа, и все члены класса, как и члены интерфейса, заключены в фигурные скобки.

public abstract class TIntegrator: IIntegrator.

{.

}.

В заголовке класса обязательно указывается служебное слово class. В данном случае указан дополнительный модификатор abstract, который означает, что класс TIntegrator будет содержать не реализованные методы. Вспомним, что у интерфейса все методы не реализованы. Абстрактный класс, как и интерфейс, не может создавать объекты. Но, в отличие от интерфейса, абстрактный класс может иметь поля и реализовывать любое количество методов и свойств. В заголовке класса TIntegrator указано так же, что он наследует от интерфейса IIntegrator. Это означает, что класс TIntegrator берет на себя обязательство реализовать все методы и свойства, описанные в интерфейсе IIntegrator. Другими словами, если наш класс TIntegrator не реализует (хотя бы в абстрактной форме) какой-либо из членов интерфейса IIntegrator, транслятор не пропустит такой код.

Рассмотрим подробнее члены класса TIntegrator.

В начале опишем ряд постоянных полей, которые потребуются для реализации некоторых методов класса TIntegrator и его наследников.

public const double minTolerance = 1e-15;

public const double maxTolerance = 0.01;

public const double defTolerance = 1e-6;

protected internal const double epsilon = 1e-16;

Постоянные поля имеют модификатор const. Их значения определяются на стадии конструирования класса как типа в том смысле, что эти значения фиксированы для всех объектов класса и его наследников, если имеют достаточный уровень доступа.

Первые два поля minTolerance, maxTolerance ограничивают допустимые значения, которые может иметь погрешность интегрирования дифференциальных уравнений. Поле defTolerance определяет погрешность по умолчанию, то есть значение погрешности, которое будет приниматься в любом объекте интегратора, если это значение не задано явно.

Эти три постоянных поля имеют максимально высокий уровень доступа public. Они видны не только методам класса TIntegrator и его наследников, но и всем приложениям, использующим библиотеку, в которой описан класс TIntegrator. В любой части кода можно написать, в частности, выражение.

double d = TIntegrator. maxTolerance;

Выполнение этого кода поместит значение постоянного поля maxTolerance в объект d.

Однако запись вида.

TIntegrator.maxTolerance = 0.1;

приведет к ошибке компиляции, так как значение постоянного поля не может быть изменено нигде и никем.

Постоянное поле epsilon имеет смысл минимальной по сравнению с единицей величины, которая используется в некоторых методах класса TIntegrator и его наследников. Поле epsilon имеет уровень доступа protected internal. Уровень доступа protected internal позволяет использовать поле epsilon только в самом классе TIntegrator и его наследниках, описанных в той же библиотеке, что класс TIntegrator.

Далее будем описывать переменные поля и свойства класса TIntegrator.

Следующим кодом класс TIntegrator реализует два свойства N и Equations интерфейса IIntegrator.

protected internal ushort n;

public ushort N { get { return n; } }.

protected internal TEquations equations;

public TEquations Equations { get { return equations; } }.

Как мы видим, поля n и equations, хранящие текущие значения числа уравнений и ссылку на метод, возвращающий уравнений, имеют модификатор доступа protected internal, более низкий, чем свойства. Причем свойства N и Equations позволяют лишь читать значения этих полей, но не менять их. Такая практика программирования препятствует без нужды рисковать значениями важных полей в период эксплуатации объектов.

Метод доступа get в том и другом свойстве содержит характерный оператор return. Оператор return возвращает значение метода-функции, в данном случае метода get, и прерывает выполнение метода.

Далее можно описать поля.

protected internal double [ ] currY;

protected internal double currt;

Они предназначены для хранения текущего (current) вектора неизвестных функций currY и текущего значения независимой переменной currt. Поле currY описано как ссылка на массив объектов типа double. На это указывают квадратные скобки. Длина массива не указана. Она устанавливается после задания числа уравнений и будет равна числу уравнений n. Нумерация элементов массива начинается индексом 0, поэтому фактически поле currY будет содержаться ссылка на массив с элементами currY[0],…, currY[n — 1].

Опишем еще два поля и соответствующие им свойства.

protected internal double stepSize;

public double StepSize { get { return stepSize; } }.

protected internal ulong equationsCallNmb;

public ulong EquationsCallNmb { get { return equationsCallNmb; } }.

Поле stepSize хранит текущий шаг интегрирования. Оно доступно только наследникам класса TIntegrator, описанным в той же библиотеке классов (модификатор доступа protected internal говорит об этом). Его свойство StepSize доступно всем (имеет модификатор доступа public). Но свойство StepSize лишь возвращает значение поля stepSize, и не может его изменить.

Поле equationsCallNmb хранит показания счетчика вызова уравнений, который может потребоваться пользователю для оценки эффективности того или иного метода интегрирования. Тип ulong определяет объекты длинных целых чисел без знака, занимающих 8 байт памяти. Пределы значений таких чисел от 0 до 18,446,744,073,709,551,615. Доступное всем свойство EquationsCallNmb возвращает пользователю значение поля equationsCallNmb.

Теперь опишем метод, возвращающий интегратор в состояние начала интегрирования.

protected internal virtual void Reset ().

{.

stepSize = Double. NaN; equationsCallNmb = 0;

}.

Модификатор virtual указывает на то, что метод Reset является виртуальным. Виртуальные методы могут быть перекрыты (изменены) в классах-наследниках. Служебное слово void (пустой тип) перед именем метода Reset означает, что метод после своего выполнения не возвращает объект какого-либо типа. Такой тип метода в C# аналогичен процедуре в Паскале. Обратите внимание и на то, что хотя метод Reset не имеет параметров, синтаксис языка C# требует написания круглых скобок после имени метода.

В теле метода находится два простых оператора присваивания. Первый оператор присваивает полю stepSize значение Double.NaN. NaN это постоянное поле класса Double (служебное слово double является псевдонимом класса Double), которое указывает на особое значение Not a Number объекта типа double. Фактически это неопределенное значение вещественного числа типа double, или представление числа «не имеющего значения». До начала интегрирования размер шага stepSize не определен, так как шаг зависит от начальных условий и самих дифференциальных уравнений.

Таким образом, метод Reset делает размер шага stepSize не определенным и обнуляет значение счетчика обращений к уравнениям equationsCallNmb.

Теперь реализуем индексатор this, описанный в интерфейсе IIntegrator.

public double this [ushort index].

{.

get.

{.

return index == 0? currt: currY [index — 1];

}.

set.

{.

if (index == 0) currt = value;

else currY [index — 1] = value;

// Любое внешнее изменение полей currt или currY.

// предполагает, что задаются начальные условия.

// Это должно приводить к вызову метода Reset,.

// «стирающего» оптимальный шаг интегрирования.

if (!Double.IsNaN (stepSize)) Reset ();

}.

}.

Выражение.

index == 0? currt: currY [index — 1].

из метода get называется условным. В нем находится оператор сравнения ==, означающий «равен ли объект index нулю»? Оператор сравнения возвращает логическое значение true или false. Если это true, то условное выражение возвращает то, что следует за знаком вопроса? (в данном случае currt). В противном случае условное выражение возвращает то, что следует за знаком двоеточия: (в данном случае элемент массива currY [index — 1]). Таким алгоритмом индексатор класса возвращает при нулевом индексе текущее значение независимой переменной currt, а при других значениях текущее значение соответствующей неизвестной функции.

На практике такое описание индексатора приведет к следующему эффекту. Пусть описан объект integrator класса TIntegrator (вернее, наследника класса TIntegrator, так как сам класс TIntegrator абстрактный и не может создавать объектов). Тогда использование индексатора может иметь вид оператора.

double d = integrator[5];

Это означает, что объекту d будет передано значение, возвращаемое методом get индексатора, отвечающее значению 5 параметра index. Согласно приведенному выше описанию индексатора, это должен быть 4-ый элемент массива неизвестных функций currY.

Но что произойдет, если у интегратора задано всего 4 уравнения и массив currY имеет границы от 0 до 3? Произойдет останов по ошибке. Возникнет, как говорят, исключительная ситуация, связанная с тем, что совершена попытка выбрать элемент массива вне зоны его границ.

Как вариант, можно изменить содержание метода get. Усилим условие, накладываемое на значение параметра index, подставив во вторую часть условного выражения еще одно условное выражение вида.

index<=n? currY [index — 1]: Double.NaN.

Полный оператор метода get теперь будет иметь вид.

return index == 0? currt: index<=n? currY [index — 1]: Double. NaN;

Теперь значение, которое будет возвращать индексатор при любом индексе, превышающем число уравнений, будет неопределенным (Double.NaN). Исключительная ситуация не возникнет.

Перейдем к методу set индексатора this.

Метод set содержит условные операторы типа if… else и типа if. В операторе if условие формулируется в круглых скобках. В данном случае в первом операторе if это условие равенства нулю индекса index == 0. Во втором операторе if условие! Double. IsNaN (stepSize) выполняется, если поле stepSize является определенным, то есть не равно NaN. Метод IsNaN класса Double возвращает логическое значение true или false в зависимости от того, имеет ли параметр этого метода значение NaN или нет. Унарный оператор отрицания ! превращает true в false и наоборот.

При выполнении условия index == 0 метод set передает полю текущего значения независимой переменной currt значение своего параметра value. В противном случае, когда индекс не равен нулю, выполняется оператор, следующий за словом else. Здесь значение параметра value метода set передается соответствующему элементу массива текущих значений неизвестных функций currY [index — 1].

В пользовательском классе метод set индексатора будет вызываться в операторе вида (это пример!).

integrator[5] = 0.1;

Вызов этого оператора фактически сводится к вызову метода set индексатора с параметром value, равным 0.1, и параметром index, равным 5.

Если индексы массива неизвестных функций currY ограничены, как и в предыдущем примере, значениями 0…3 (система содержит четыре уравнения), то вызов метода set приводит к возникновению исключительной ситуации.

Для того, чтобы избежать исключительной ситуации, можно усилить условие, дополнительно потребовав, чтобы значение индекса не превышало число уравнений n. Это усложнит условный оператор if… else, придав ему вид.

if (index == 0) currt = value;

else if (index <= n) currY [index — 1] = value;

Здесь, как видно, при условии index > n метод set вообще не выполняет никаких действий. Поэтому в новом варианте кода оператор integrator[5] = 0.1 не выполнит никаких действий, если число уравнений меньше пяти.

Тем и хороши свойства, что они могут оберегать поля от неразумного использования.

Заметим, правда, что в описанных реализациях методов set и get нельзя с уверенностью сказать, какой из вариантов кода предпочтительнее. Если пользователь получит сообщение об исключительной ситуации, то он разберется, поймет ошибку и исправит ее. Во втором варианте код выполнится без помех и пользователь будет введен в заблуждение, считая, что его код не содержит ошибочных операторов.

В последнем операторе if метода set вызывается метод Reset, если поле stepSize имеет осмысленное значение (не NaN). Отличие stepSize от NaN фактически означает, что метод Reset еще не вызывался. Тем самым метод Reset вызывается лишь один раз при попытке изменить значение независимой переменной, либо одной из неизвестных функций. При любом таком изменении интегратор предполагает, что пользователь этим изменением задает новые начальные условия интегрирования.

Продолжим рассмотрение других членов класса TIntegrator.

protected internal double tolerance;

public virtual double Tolerance.

{.

get { return tolerance; }.

set.

{.

tolerance =.

value > maxTolerance? maxTolerance :

value < minTolerance? minTolerance :

value;

}.

}.

Так описывается поле tolerance, хранящее текущую погрешность интегрирования, и соответствующее ему свойство Tolerance.

Свойство Tolerance скрытого поля tolerance доступно всем. Оно реализует абстрактное описание свойства Tolerance интерфейса IIntegrator. Модификатор virtual в описании свойства означает, что оно является виртуальным и его методы set и get могут быть перекрыты, изменены в классах — наследниках класса TIntegrator.

Обратим внимание на содержание метода set свойства Tolerance. Если пользователь задает погрешность, выходящую за границы, определенные постоянными полями maxTolerance и minTolerance, то поле tolerance принимает значение той границы, за которую выходит значение параметра value.

Возможен другой вариант кода, предусматривающий в этом случае возбуждение исключительной ситуации. Тогда пользователь предупреждается о нарушении границ погрешности. В этом варианте тело метода set может выглядеть следующим образом.

if (value > maxTolerance || value < minTolerance).

{.

throw (new ApplicationException (.

" Погрешность выходит за допустимые границы."));

}.

tolerance = value;

Теперь выполняется полная проверка значения параметра value, с которым пользователь передает предполагаемую погрешность. Если это значение больше максимально допустимого maxTolerance или (оператор || означает «или» в C#) меньше минимально допустимого minTolerance, то стандартный оператор throw создает («вбрасывает») исключительную ситуацию. Параметром оператора throw является объект стандартного класса ApplicationException. Этот объект создается вызовом оператора new и конструктора класса ApplicationException с параметром-строкой «Погрешность выходит за допустимые границы.». Эта строка будет выведена как комментарий, или сообщение об ошибке, отвечающей возникшей исключительной ситуации. После оператора throw выполнение метода set прерывается, и программа ищет обработчик исключительной ситуации, готовый принять созданный объект с его сообщением. Обработчик устанавливается пользователем в расчете на возможную исключительную ситуацию. С технологией обработки исключительных ситуаций мы познакомимся несколько позже.

Опишем еще одно поле и еще один метод класса TIntegrator.

protected internal double [ ] startLocalVector;

protected internal virtual void Initialize ().

{.

// Создается экземпляр объекта currY и.

// резервируется память под массив текущих значений неизвестных функций.

currY = new double [n];

// Элементы массива currY делаются неопределенными (NaN, Not a Number).

for (ushort i = 0; i < n; i++) currY [i] = Double. NaN;

// Независимая переменная делается неопределенной NaN.

currt = Double. NaN;

// Создается экземпляр объекта startLocalVector и.

// резервируется память под массив правых частей уравнений.

startLocalVector = new double [n];

// Погрешность принимает значение «по умолчанию».

Tolerance = defTolerance;

}.

Поле-массив startLocalVector используется в некоторых методах класса TIntegrator и его наследников. Этот массив хранит начальное значение вектора правой части дифференциальных уравнений, необходимый для оценки начального шага интегрирования.

В виртуальном методе Initialize с ограниченным доступом объявленные выше поля currY, currt, startLocalVector, и tolerance принимают начальные значения. Предполагается, что к моменту вызова метода Initialize число уравнений n установлено. Поля-ссылки на массивы currY и startLocalVector инициализируются с помощью оператора new и конструктора double [n]. Оператор выделяет память под хранение необходимого числа элементов массива и передает ссылку на эту память указанному полю. Конструктор сообщает всю необходимую информацию о создаваемом объекте оператору new.

Обратите внимание на оператор цикла for (;;). Вслед за ним находится один или несколько операторов, объединенных фигурными скобками. В нашем коде это один оператор присваивания currY [i] = Double.NaN. Этот оператор называют «оператором, которым управляет цикл», или, короче, «телом цикла».

Синтаксис языка C# требует, чтобы оператор цикла типа for (;;) имел в круглых скобках три выражения, разделенных точкой с запятой. В первом выражении параметр цикла принимает начальное значение и может быть описан (как в нашем примере). Второе выражение является условием выполнения тела цикла (в нашем случае условие выполнения тела i < n). Наконец, в третьем выражении параметр цикла изменяется. В нашем операторе изменение параметра определяется выражением i++. В языке C# операция i++ означает «добавление к i единицы и сохранение нового значения в объекте i».

Оператор for (;;) работает следующим образом:

сразу после установки начального значения проверяется условие выполнения тела. Если оно соблюдается, то тело выполняется, если нет, то тело цикла не выполняется, а совершается «выход из цикла», то есть переход к оператору, следующему за оператором цикла вместе с его телом.

После выполнения тела цикла управление передается разделу, изменяющему параметр цикла. Параметр цикла меняется.

Затем вновь проверяется условие выполнения тела цикла, и т. д.

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

Продолжим описание членов класса TIntegrator.

protected internal sbyte direction;

public sbyte Direction { get { return direction; } }.

protected internal CallBackEventHandler callBack;

public virtual CallBackEventHandler CallBack.

{ get { return callBack; } set { callBack = value; } }.

Поле direction хранит направление изменения независимой переменной (плюс один при росте и минус один при уменьшении). Соответствующее свойство Direction доступно всем. Оно возвращает значение поля direction. Тип sbyte означает, что объект этого типа хранит целое число со знаком объемом в один байт. Пределы значений объектов типа sbyte [-128; 127].

Поле callBack хранит ссылку на делегата типа CallBackEventHandler, передающего от пользователя обработчики, или функции обратного вызова. Об этом шла речь в описании интерфейса IIntegrator. Свойство CallBack реализует свойство интерфейса IIntegrator с тем же именем и является, как мы видим, виртуальным. Наследник может изменить содержание его методов set и get.

Далее разместим и прокомментируем описание конструкторов объектов класса TIntegrator.

protected internal TIntegrator (.

ushort N,.

TEquations Equations,.

CallBackEventHandler CallBack.

).

{.

if (N == 0).

throw (new ApplicationException («Число уравнений равно нулю!?»));

n = N;

if (Equations == null).

throw (new ApplicationException («Уравнения не заданы!!!»));

equations = Equations;

Initialize ();

this.CallBack = CallBack;

}.

protected internal TIntegrator (ushort N, //Число уравнений.

TEquations Equations//уравнения.

): this (N, Equations, null) { }.

В языке C# конструкторы имеют те же имена, что и имя класса. По виду это обычные методы с заголовками, содержащими модификатор доступа, имя и список параметров в круглых скобках. В данном случае модификатор доступа описанных конструкторов ограничен наследниками, описанными в той же библиотеке (protected internal), так как класс TIntegrator абстрактный и не может создавать объекты. Ему не нужны доступные всем конструкторы. Важное отличие описания конструктора от обычного метода в том, что конструктор не возвращает никакой тип, даже void. Описание обычного метода должно указывать на тип, возвращаемого объекта.

В классе TIntegrator описаны два конструктора. Но число параметров у них разное. Такой способ описания нескольких методов с одинаковыми именами, но разной сигнатурой, называется перегрузкой (overloading).

У первого конструктора три параметра N, Equations и CallBack. Они определяют соответственно число уравнений, делегат уравнений и делегат функции обратного вызова. После проверки значений этих параметров условными операторами конструктор присваивает их соответствующим полям. Если число уравнений равно нулю, или уравнения не заданы (это означает, что ссылка Equations равна null — это ссылка в «никуда»), то возбуждаются исключительные ситуации. В строке сообщения при создании объекта класса исключительной ситуации ApplicationException присутствует символ. Это пример так называемой escape-последовательности, используемой для вывода специальных символов. Начинается escape-последовательность с символа (back slash). За этим следуют другие буквы из некоторого ограниченного набора. Так символ n означает, что перед выводом следующей далее строки в окно необходимо перейти на новую строку (new line).

Далее вызывается метод Initialize, рассмотренный выше. Поле callBack инициализируется через свое виртуальное свойство CallBack, а не напрямую.

Второй конструктор имеет только два параметра. С его помощью можно создавать объект, которому не нужна ссылка на функцию обратного вызова CallBack. Второй конструктор оператором this (N, Equations, null), отделенным двоеточием, просто вызывает первый конструктор, подставляя в качестве третьего параметра значение null — обработчик CallBack отсутствует. Тело второго конструктора пустое, но скобки {} писать необходимо.

Рассмотрим следующий метод класса TIntegrator.

protected internal virtual bool DoCallBack ().

{.

CallBackEventArgs e = new CallBackEventArgs ();

callBack (this, e);

return e. Stop;

}.

Виртуальный метод DoCallBack возвращает объект типа bool. Так обозначается логический тип, принимающий два значения true и false. Метод DoCallBack используется классом TIntegrator для вызова обработчика, или функции обратного вызова callBack, заданной пользователем.

В теле метода DoCallBack создается экземпляр e класса CallBackEventArgs, который должен содержать состояние события, обрабатываемого пользователем. Фактически у объекта e есть единственное поле Stop логического типа. Поле Stop может быть изменено пользователем, после того как будет вызван метод callBack.

При вызове метода callBack первым параметром является this. Это ссылка на объект класса, в котором произошло событие. Таковым является класс TIntegrator или его наследник. Метод DoCallBack должен вызываться классом TIntegrator в точке события.

После вызова обработчика callBack в методе DoCallBack стоит оператор возврата return, который прекращает выполнение метода DoCallBack и одновременно возвращает значение поля e. Stop как результат функции DoCallBack. Таким образом, если метод DoCallBack возвращает true, то цикл интегрирования должен прервать свое выполнение, а если false, то цикл интегрирования должен продолжиться.

Далее в классе находятся описания трех методов, реализующих поиск оптимального начального шага интегрирования. Они содержат алгоритм такого поиска, в детали которого мы не будем вдаваться. Их полный код приведен в тексте файла с описанием класса TIntegrator. Прокомментируем лишь заголовки и некоторые детали использования этих методов.

protected internal virtual void FirstEstimate (ref double df, ref double dy).

{.

…

}.

protected internal virtual void SecondEstimate.

(ref double der2, double [ ] tempLocalVector).

{.

…

}.

double StartStepSize (double MaxStepSize).

{.

double df = 0.0, dy = 0.0, initStep;

FirstEstimate (ref df, ref dy);

…

double der2 = 0.0, der12, h1;

SecondEstimate (ref der2, tempLocalVector);

…

}.

Методы FirstEstimate и SecondEstimate являются виртуальными. Они обеспечивают алгоритмы первичной и вторичной оценки (estimation) шага интегрирования и вызываются внутри метода StartStepSize.

Обратим внимание на модификатор ref, использованный при описании параметров методов FirstEstimate и SecondEstimate. Модификатор ref означает, что вызов объекта (в данном случае типа double) производится по ссылке (reference). По правилу языка C# фактический параметр, который подставляется в метод при его вызове на место параметра, вызываемого по ссылке, должен так же иметь модификатор ref. Именно это указано в операторах тела метода StartStepSize, вызывающего методы FirstEstimate и SecondEstimate. При этом фактический параметр с модификатором ref должен к моменту обращения иметь определенное значение. В данном случае параметры df и dy внутри метода StartStepSize принимают нулевые значения перед вызовом метода FirstEstimate. Так же нулевое значение принимает параметр der2 перед вызовом SecondEstimate. Значения параметров по ссылке могут быть изменены вызываемым методом. Другими словами значения df и dy могут отличаться от нуля после вызова метода FirstEstimate. Так же вызов метода SecondEstimate может изменить значение объекта der2.

Метод StartStepSize возвращает объект типа double. Это будет найденный им оптимальный шаг интегрирования. Метод StartStepSize имеет один параметр MaxStepSize, который используется при определении оптимального шага. Это параметр по значению и, в отличие от параметров по ссылке, значение MaxStepSize метод StartStepSize изменить не может.

В описании метода StartStepSize отсутствует модификатор доступа. Это означает, что модификатором является private. Если модификатор доступа не указывается при описании какого-либо члена класса (описание «по умолчанию»), то им является модификатор private. Это правило языка. Член класса с модификатором private доступен только членам этого же класса, и никому более.

Прокомментируем два оставшихся метода класса.

public bool IntegrateTo (double tEnd).

{.

for (ushort i = 0; i < n + 1; i++).

if (Double.IsNaN (this [i])).

throw (new ApplicationException.

(«Не заданы начальные условия интегрирования!»));

if (tEnd — currt == 0).

throw (new ApplicationException («Шаг интегрирования равен нулю?!»));

direction = (sbyte)(tEnd — currt > 0? 1: -1);

try.

{.

Equations (currt, currY, startLocalVector);

equationsCallNmb++;

}.

catch (Exception e).

{.

throw (new ApplicationException.

(e.Message + «Ошибка при первом обращении к уравнениям!»));

}.

if (Double.IsNaN (stepSize)).

{.

stepSize = StartStepSize (Math.Abs (tEnd — currt));

if (Double.IsNaN (stepSize) || Double. IsInfinity (stepSize)).

{.

throw (new ApplicationException.

(«Ошибка при вычислении начального шага интегрирования!»));

}.

}.

else stepSize = Math. Abs (stepSize) * Direction;

return BasicLoop (tEnd);

}.

protected internal abstract bool BasicLoop (double tEnd);

Метод IntegrateTo реализует абстрактный метод интерфейса IIntegrator с тем же именем. У метода IntegrateTo есть один параметр по значению tEnd, в который подставляется значение независимой переменной, до которого следует интегрировать систему уравнений. Метод IntegrateTo возвращает объект типа bool, равный true, если метод завершен «естественным путем», то есть конечная точка интегрирования tEnd достигнута. В самом конце тела метода IntegrateTo стоит оператор return BasicLoop (tEnd), вызывающий метод BasicLoop и возвращающий то значение, которое возвращает BasicLoop.

Описание метода BasicLoop содержит модификатор abstract. Поэтому метод BasicLoop не имеет тела и, собственно поэтому, класс TIntegrator является абстрактным. Метод BasicLoop должен быть реализован в наследниках класса TIntegrator. Там он будет содержать конкретный алгоритм интегрирования системы дифференциальных уравнений.

Рассмотрим последовательно тело метода IntegrateTo.

При создании объекта класса TIntegrator и его наследников вызывается конструктор. В свою очередь конструктор вызывает метод Initialize, который придает неопределенные значения (NaN) независимой переменной currt и неизвестным функциям currY. Эти значения сохраняются в полях объекта до тех пор, пока их не изменит обращение к индексатору this. Своим методом set индексатор изменит значения currt и currY на те числа, которые будут введены в качестве начальных условий. Предполагается, что только после этого объект будет вызывать метод интегрирования IntegrateTo.

Первые операторы метода IntegrateTo проверяют, заданы ли начальные условия. Для этого просто проверяется, определены ли поля currt и currY. При этом вызывается метод IsNaN класса Double, параметром которому передается индексатор this в цикле при всех допустимых значениях индекса от 0 до n. Если хотя бы одно из проверяемых значений не определено, метод IsNaN возвратит true. Это, как указано в коде, приведет к возбуждению исключительной ситуации типа ApplicationException с соответствующим сообщением и цикл прервет свое выполнение, как, впрочем, и весь метод IntegrateTo. Код передаст созданный объект исключительной ситуации первому попавшемуся обработчику.

Если начальные условия заданы, то прерывание не наступает. Далее начальное значение независимой переменной currt сравнивается с ее конечным значением tEnd — параметром метода IntegrateTo. При совпадении этих значений так же создается объект исключительной ситуации с сообщением «Шаг интегрирования равен нулю?!» и выполнение метода IntegrateTo прерывается.

Следующий оператор

direction = (sbyte)(tEnd — currt > 0? 1: -1);

устанавливает значение поля direction — направление изменения независимой переменной при интегрировании. Обратите внимание на преобразование типов целых переменных. Объекты 1 или -1 имеют тип int. Но поле direction имеет тип sbyte. Для перехода от объектов типа int к объектам типа sbyte необходимо писать явный оператор преобразования типов (sbyte), как указано в коде.

Рассмотрим подробнее следующий оператор вида.

try.

{.

Equations (currt, currY, startLocalVector);

equationsCallNmb++;

}.

catch (Exception e).

{.

throw (new ApplicationException.

(e.Message + «Ошибка при первом обращении к уравнениям!»));

}.

Это пример оператора, в котором производится обработка исключительной ситуации.

Блок, выделенный служебным словом try (пытайся) содержит код, который должен быть выполнен. В данном случае здесь вызывается метод Equations, вычисляющий вектор правой части дифференциальных уравнений startLocalVector. в точке currt, currY. Затем значение счетчика обращения к уравнениям equationsCallNmb увеличивается на единицу. Если при выполнении этих операций возникнет исключительная ситуация (скорее всего она может возникнуть только внутри метода Equations, заданном пользователем), то выполнение дальнейших операторов будет прервано и управление будет немедленно передано в блок catch.

В блоке catch происходит обработка исключительной ситуации. Объект исключительной ситуации попадает в блок catch через параметр e типа Exception. Это стандартный класс из библиотеки .NET, который является предком всех классов, порождающих объекты исключительных ситуаций. В частности, уже встречавшийся ранее класс ApplicationException является наследником Exception. Каждый объект класса Exception и его наследников содержит поле Message, несущее сообщение об исключительной ситуации. В блок catch объект прибывает с конкретным значением поля Message. В нашем коде, при обработке исключительной ситуации внутри блока catch мы создаем новый объект исключительной ситуации типа ApplicationException. В качестве сообщения мы переправляем новому объекту сообщение e. Message о том, что произошло в методе Equations, и добавляем свое сообщение «Ошибка при первом обращении к уравнениям!» .

Следующий условный оператор в методе IntegrateTo.

if (Double.IsNaN (stepSize)).

проверяет, имеет ли текущий шаг интегрирования stepSize осмысленное значение. Значение NaN поля stepSize означает, что шаг еще не определен. Тогда вызывается метод StartStepSize.

stepSize = StartStepSize (Math.Abs (tEnd — currt));

Он вычисляет оптимальный начальный шаг интегрирования. Параметром метода StartStepSize является максимальный шаг — длина всего интервала интегрирования от начальной точки currt до конечной tEnd.

В особой ситуации метод может возвратить неопределенное StartStepSize значение шага. Полученный шаг может иметь значение NaN или быть бесконечным. Если это происходит, то следующий условный оператор возбуждает исключительную ситуацию, чтобы прервать выполнение метода IntegrateTo.

if (Double.IsNaN (stepSize) || Double. IsInfinity (stepSize)).

{.

throw (new ApplicationException.

(«Ошибка при вычислении начального шага интегрирования!»));

}.

Здесь метод IsInfinity класса Double проверяет, имеет ли аргумент (шаг stepSize) бесконечное значение.

Последний оператор метода IntegrateTo вызывает метод BasicLoop, который должен содержать алгоритм интегрирования системы дифференциальных уравнений. В классе TIntegrator метод BasicLoop описан как абстрактный и не имеет тела.

protected internal abstract bool BasicLoop (double tEnd);

Это последний член класса TIntegrator. На этом описание класса TIntegrator завершается. Можете посмотреть весь текст модуля с описанием класса TIntegrator. Обратите внимание на структуру комментариев и операторы, из которых состоят тела не рассмотренных методов FirstEstimate, SecondEstimate и StartStepSize. В теле этих методов встречаются обращение к стандартным математическим функциям, сосредоточенным в классе Math библиотеки .NET. Это такие функции, как абсолютное значение Abs, квадратный корень Sqrt, степень числа Pow, минимальное Min и максимальное Max из двух чисел. Алгоритм вычисления начального оптимального шага взят из работы E. Hairer, S.P. Norsett and G. Wanner, Solving ordinary differential equations I, nonstiff problems, 2nd edition, Springer Series in Computational Mathematics, Springer-Verlag (1993).

http://scitation.aip.org/getabs/servlet/GetabsServlet?prog=normal&id=SIREAD000032000003000485000001&idtype=cvips&gifs=yes.

Класс RK4.

В качестве следующего примера рассмотри класс интегратора, построенный на основе алгоритма Рунге-Кутта 4-ого порядка и позволяющий создавать объекты. Назовем новый класс RK4. Наследуем RK4 от абстрактного класса TIntegrator.

В начале приведем общий список членов класса RK4 с комментариями В заголовке класса RK4 через двоеточие указано, что он наследует от класса TIntegrator. Следовательно, все члены класса TIntegrator являются так же членами его наследника RK4. В описании класса RK4 указываются лишь новые члены и те методы и свойства, которые класс RK4 перекрывает.

public class RK4: TIntegrator.

{.

/*Класс RK4 добавляет новые поля, описанные без модификатора доступа и, следовательно, имеющие доступ private — только для членов класса RK4. Это одно поле lowestAcc типа double и пять рабочих массивов tempY, tempLocalVector, oneStepLocalVector, firstHalfStepLocalVector, secondHalfStepLocalVector, элементами которых являются также объекты типа double.

Свойство lowestAcc хранит пороговое значение погрешности. Это свойство используется в алгоритме выбора оптимального шага. Во время интегрирования шаг может оказаться излишне мелким. Таким, что погрешность станет ниже порогового значения. При этом условии алгоритм увеличит шаг для оптимизации счета. Поля-массивы используются для организации алгоритма счета в методах класса.*/.

double lowestAcc;

double [ ] tempY, tempLocalVector, oneStepLocalVector,.

firstHalfStepLocalVector, secondHalfStepLocalVector;

/*Ниже описано виртуальное свойство погрешности Tolerance. В описании свойства Tolerance присутствует модификатор доступа override. Этим указывается, что виртуальное свойство Tolerance перекрыто в классе RK4. Внутри методов set и get свойства Tolerance есть оператор, вызывающий методы set и get предка. Это base.Tolerance. Как видно метод get свойства Tolerance просто вызывает метод get предка. В перекрытом методе set проводится дополнительная установка значения пороговой погрешности lowestAcc. */.

public override double Tolerance.

{.

get { return base. Tolerance; }.

set.

{.

base.Tolerance = value;

// Пороговая погрешность.

lowestAcc = 0.025 * tolerance;

}.

}.

/*Дальше описывается метод Initialize, который так же перекрывает виртуальный метод Initialize класса TIntegrator. Код метода Initialize рассмотрим несколько ниже.*/.

protected internal override void Initialize ().

{.

//…

}.

/*Конструкторы класса-предка должны быть обязательно перекрыты в наследниках, даже, если они не вносят никакого дополнительного кода. Это правило языка. Обратите внимание, что конструкторы класса RK4 описаны с модификатором public. Класс RK4 не является абстрактным — он не содержит не реализованных свойств и методов. Поэтому его конструкторы должны быть доступны.

В описании конструкторов RK4 следует обратить внимание на операторы вида base (N, Equations), отделенные двоеточием от имени конструктора. Это вызов конструкторов предка. Тела новых конструкторов RK4 в данном случае пусты, но фигурные скобки писать необходимо. */.

public RK4 (ushort N, TEquations Equations): base (N, Equations) { }.

public RK4 (ushort N, TEquations Equations,.

CallBackEventHandler CallBack).

: base (N, Equations, CallBack) { }.

/*Класс вводит новый метод DoStep, который реализует формулу Рунге-Кутта 4-ого порядка. Тело метода DoStep будет прокомментировано ниже. Здесь лишь отметим, что этот метод доступен только классам-наследникам из той же библиотеки (модификатор доступа protected internal) и он является виртуальным.*/.

protected internal virtual void DoStep (double tStart, double currStep,.

double [ ] yStart, double [ ] yEnd).

{.

// Метод Рунге-Кутта 4-ого порядка.

}.

// Класс реализует метод BasicLoop, содержащий цикл интегрирования.

protected internal override bool BasicLoop (double tEnd).

{.

//…

}.

}.

Метод Initialize, реализованный классом RK4, добавляет к унаследованной версии инициализацию рабочих массивов-полей.

protected internal override void Initialize ().

{.

base.Initialize ();

try.

{.

oneStepLocalVector = new double [n];

firstHalfStepLocalVector = new double [n];

secondHalfStepLocalVector = new double [n];

tempY = new double [n];

tempLocalVector = new double [n];

}.

catch (Exception e).

{.

throw (new ApplicationException.

(e.Message + «Прерывание в методе Initialize.»));

}.

}.

Первый оператор метода base. Initialize (); вызывает унаследованную версию метода Initialize. Обратите внимание на блоки try и catch, предназначенные для обнаружения и обработки исключительной ситуации в том случае, если рабочие массивы займут слишком много памяти.

Метод DoStep реализует формулу Рунге-Кутта 4-ого порядка (см. Г. Корн, Т. Корн, Справочник по математике, стр. 702, Таблица 20.8−1, формула c), которая имеет вид.

yEnd = yStart + (k½ + k2+ k3+ k4/2) / 3, где.

k½ = f (tStart, yStart)*h/2;

k2 = f (tStart + h/2, yStart + k½)*h;

k3 = f (tStart + h/2, yStart + k2/2)*h;

k4/2 = f (tStart + h, yStart + k3)*h/2;

Здесь функция f является вектором правой части уравнений yi = fi (t, y1,…, yN). Ссылка на этот вектор возвращается методом equations в третьем параметре (tempLocalVetor в методе DoStep). Независимая переменная tStart, вектор неизвестных функций yStart, шаг h, обозначенный currStep, и возвращаемый вектор yEnd являются параметрами метода DoStep. Сам метод DoStep вычисляет вектор yEnd, последовательно применяя приведенные формулы.

protected internal virtual void DoStep (double tStart, double currStep,.

double [ ] yStart, double [ ] yEnd).

{.

// Метод Рунге-Кутта 4-ого порядка.

// Реализуется формула yEnd=yStart+(k½+k2+k3+k4/2)/3.

double HalfStep = 0.5 * currStep;

try.

{.

// tempLocalVector=f (tStart, yStart);

// k½ =tempLocalVector*currStep/2.

equations (tStart, yStart, tempLocalVector);

equationsCallNmb++;

for (ushort count = 0; count < n; count++).

{.

// k½.

yEnd [count] = HalfStep * tempLocalVector [count];

// yStart + k½.

tempY [count] = yStart [count] + yEnd [count];

};

// tempLocalVector=f (tStart+currStep/2,yStart+k½);

// k2=tempLocalVector*currStep.

equations (tStart + HalfStep, tempY, tempLocalVector);

equationsCallNmb++;

for (ushort count = 0; count < n; count++).

{.

// yS+k2/2.

tempY [count] = yStart [count] + HalfStep * tempLocalVector [count];

// k½+k2.

yEnd [count] = yEnd [count] + currStep * tempLocalVector [count];

};

// tempLocalVector=f (tStart+currStep/2,yStart+k2/2);

// k3=tempLocalVector*currStep.

equations (tStart + HalfStep, tempY, tempLocalVector);

equationsCallNmb++;

for (ushort count = 0; count < n; count++).

{.

//yStart+k3.

tempY [count] = yStart [count] + currStep * tempLocalVector [count];

//k½+k2+k3.

yEnd [count] = yEnd [count] + currStep * tempLocalVector [count];

};

// tempLocalVector=f (tStart+currStep, yStart+k3);

// k4/2=tempLocalVector*currStep/2.

equations (tStart + currStep, tempY, tempLocalVector);

equationsCallNmb++;

}.

catch (Exception e).

{.

throw (new ApplicationException.

(e.Message + «Ошибка при обращении к уравнениям!»));

}.

//Результирующие значения функций после данного шага.

for (ushort count = 0; count < n; count++).

{.

// yEnd = yStart+(k½+k2+k3+k4/2)/3.

yEnd [count] = yStart [count] + 0.333 333 333 333 333 312 *.

(yEnd [count] + HalfStep * tempLocalVector [count]);

if (Double.IsNaN (yEnd [count]) || Double. IsInfinity (yEnd [count])).

throw (new ApplicationException.

(«Компонента NaN или Infinity в методе DoStep!»));

};

}.

Комментарии в тексте расшифровывают эту формулу и алгоритм ее реализации.

Обратите внимание на то, что практически все тело метода DoStep помещено в блок try… catch. Если при выполнении метода возникает исключительная ситуация (это, скорее всего, происходит во время одного из обращений к методу equations), обработчик блока catch возбуждает новую исключительную ситуацию, добавляя к сообщению e. Message о возникшей ситуации свое сообщение.

Теперь рассмотрим подробнее реализацию метода BasicLoop.

protected internal override bool BasicLoop (double tEnd).

{.

double currStep, currHalfStep, Divergence, maxDivergence;

bool isGoodAcc, isLastStep, isStepTooSmall;

try.

{.

// Текущий шаг currStep выбирается равным оптимальному.

currStep = stepSize;

// Перед входом в цикл вызывается обработчик.

if (callBack ≠ null && DoCallBack ()) return false;

// Цикл прохождения всего полного шага интегрирования.

do.

{.

// Если текущий шаг составляет более 80% всего интервала,.

// он выбирается равным интервалу tEnd-FCurrt.

// Флаг isLastStep регулирует выход из основного цикла.

if ((currt + 1.25 * currStep — tEnd) * direction >= 0).

{.

// Шаг является последним.

currStep = tEnd — currt; isLastStep = true;

}.

else isLastStep = false; //Это не последний шаг.

// Флаг isGoodAcc регулирует возможность увеличения шага,.

// если погрешность излишне мала.

isGoodAcc = true; //Допустим рост шага.

//Делается один текущий шаг.

DoStep (currt, currStep, currY, oneStepLocalVector);

//Цикл выбора шага, удовлетворяющего требуемой точности.

do.

{ //Делаются две половины шага.

currHalfStep = 0.5 * currStep;

DoStep (currt, currHalfStep, currY, firstHalfStepLocalVector);

DoStep (currt + currHalfStep, currHalfStep,.

firstHalfStepLocalVector, secondHalfStepLocalVector);

// Вычисляется макс. относительное расхождение.

// в результатах одного шага и двух половинок.

maxDivergence = 0.0;

for (ushort count = 0; count < n; count++).

{.

if (Double.IsInfinity (1.0 / oneStepLocalVector [count]) &&.

Double.IsInfinity (1.0 / secondHalfStepLocalVector [count])).

continue;

if (oneStepLocalVector [count] ≠ secondHalfStepLocalVector [count]).

{.

// Абсолютное расхождение.

Divergence = Math. Abs (oneStepLocalVector [count] ;

secondHalfStepLocalVector [count]);

// Относительное расхождение вычисляется.

// при достаточно большом абсолютном.

if (Divergence > epsilon).

Divergence =.

2.0 * Divergence /.

(Math.Abs (oneStepLocalVector [count]) +.

Math.Abs (secondHalfStepLocalVector [count]));

if (Double.IsInfinity (Divergence) || Double. IsNaN (Divergence)).

throw (new ApplicationException.

(«Проблема при вычислении maxDiv»));

// Корректируется значение максимального расхождения.

if (Divergence > maxDivergence) maxDivergence = Divergence;

}.

}.

//Делаются выводы по сравнению с требуемой точностью.

if (maxDivergence > tolerance).

{.

isGoodAcc = false; // Рост шага запрещен.

currStep = currHalfStep; // Шаг уменьшается вдвое.

stepSize = currHalfStep;

isLastStep = false;

firstHalfStepLocalVector.CopyTo (oneStepLocalVector, 0);

}.

isStepTooSmall = currt == 0.0 && Math. Abs (currStep) <= epsilon.

|| Math. Abs (currStep) <= Math. Abs (currt) * epsilon;

}.

// Процесс проверки и уменьшения шага продолжается.

// до достижения требуемой точности.

while.

(maxDivergence > tolerance && !isStepTooSmall);

if (isStepTooSmall).

{ throw (new ApplicationException («Шаг слишком мал»)); }.

// Точность достигнута. Устанавливаются скорректированные значения.

// функций в конце сделанного шага.

for (ushort count = 0; count < n; count++).

currY [count] =.

secondHalfStepLocalVector [count] + 0.666 666 666 666 666 624 *.

(secondHalfStepLocalVector [count] ;

oneStepLocalVector [count]);

//Если шаг не последний.

if (!isLastStep).

{.

currt += currStep; //Изменение переменной.

//Вызывается обработчик.

if (callBack ≠ null && DoCallBack ()) return false;

if (maxDivergence < lowestAcc && isGoodAcc).

{.

//Слишком мелкий шаг удваивается.

currStep += currStep;

if (Math.Abs (currStep) > Math. Abs (tEnd — currt)).

currStep = tEnd — currt;

stepSize = currStep;

}.

}.

//Прерывание цикла, если сделан последний шаг.

else break;

}.

while (true);

//Текущее значение переменной совпадает с конечным.

currt += currStep; //tEnd.

//Вызывается обработчик.

return !(callBack ≠ null && DoCallBack ());

}.

catch (Exception e).

{.

throw (new ApplicationException.

(e.Message+ «BasicLoop прерван»));

}.

}.

В начале метод BasicLoop копирует поле текущего оптимального шага stepSize в локальную переменную currStep. Значение поля stepSize либо было инициализировано унаследованным методом IntegrateTo, либо сохранилось от предыдущего вызова того же метода.

Затем вызывается обработчик callBack, который может быть задан пользователем класса. В условном операторе логическая формула callBack ≠ null && DoCallBack () возвращает true лишь в том случае, если обработчик callBack задан (ссылка на него не равна null — ссылке в «никуда») и метод DoCallBack, вызывающий обработчик, возвращает значение true. В этом случае выполнение метода BasicLoop прерывается, а сам метод возвращает false. Если обработчик callBack не определен (тогда ссылка на него имеет значение null), либо он определен, но DoCallBack возвращает после своего выполнения значение false, выполнение метода BasicLoop продолжается.

Далее в методе BasicLoop находится оператор цикла do… while (true). Он охватывает большую часть операторов метода BasicLoop. Такой оператор (условие while (true) всегда выполняется) является бесконечным циклом, и выход из него производится стандартной командой break, прерывающей выполнение любого цикла в C#. Внутри оператора do… while происходят оценки, позволяющие провести интегрирование уравнений от значения currt независимой переменной до tEnd — параметра метода BasicLoop.

В начале цикла длина всего интервала интегрирования |tEnd — currt| сравнивается с величиной текущего шага |currStep|. Если шаг превышает 80% длины интервала, то он выбирается равным всему интервалу и флаг isLastStep, указывающий на то, что это последний шаг, устанавливается в true. В противном случае флагу isLastStep присваивается значение false.

Следующий оператор устанавливает локальный флаг (объект типа bool) isGoodAcc, регулирующий возможность увеличения шага, в true.

Затем, с помощью метода DoStep вычисляется результат оптимального шага currStep по формуле Рунге-Кутта. Этот результат возвращается в виде вектора приближенного решения oneStepLocalVector в точке currt + currStep.

Следующим оператором метода BasicLoop является новый цикл do… while, который «настраивает» шаг currStep в соответствии с допустимой погрешностью.

Для этого с помощью метода DoStep вычисляется результат шага той же длины currStep, но выполненного двумя половинными шагами. Результирующие вектора запоминаются в массивах firstHalfStepLocalVector и secondHalfStepLocalVector. Между компонентами вектора secondHalfStepLocalVector и компонентами вектора oneStepLocalVector, полученного одним полным шагом, устанавливается максимальное относительное расхождение maxDivergence. Если полученное расхождение превышает допустимую погрешность tolerance, то текущий шаг currStep (как и поле оптимального шага stepSize) уменьшается вдвое, флаги «хорошего приближения» isGoodAcc и «последнего шага» isLastStep «опускаются» (принимают значение false), а компоненты вектора последней точки искомой функции oneStepLocalVector принимают значения компонент вектора, полученного после половинного шага, firstHalfStepLocalVector.

При вычислении расхождений Divergence и maxDivergence проводится ряд оценок, позволяющих возбудить исключительную ситуацию, или обойти ее при некоторых значениях переменных. В частности стандартный оператор continue позволяет в определенной ситуации обойти все проверки, попав сразу на следующий шаг цикла.

Так как уменьшение шага не может продолжаться до бесконечности, устанавливается значение флага isStepTooSmall в зависимости от значения текущего шага.

Условием окончания цикла уменьшения шага является либо условие достижения требуемой точности maxDivergence < tolerance, либо недопустимо малое значение шага (флаг isStepTooSmall поднят).

Если шаг слишком мал, то возбуждается исключительная ситуация с необходимым сообщением.

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

Далее стоит условный оператор, в результате которого цикл, а за ним и сам метод BasicLoop будет завершен, если сделанный шаг последний на заданном интервале. Если это не так, то выполняется ряд операторов, в которых.

устанавливается новое значение поля независимой переменной currt,.

вызывается обработчик callBack (если он задан).

и при условии, что флаг isGoodAcc поднят (т.е. рост шага допустим) и максимальное расхождение maxDivergence меньше нижней границы погрешности lowestAcc, текущий шаг currStep удваивается и передает свое значение полю оптимального шага stepSize. При этом если удвоенный шаг выводит точку за границу интервала tEnd, значение currStep корректируется.

Далее цикл возвращается к началу, имея новые значения currt, соответствующий вектор currY и, возможно, новое значение текущего шага currStep и поля оптимального шага stepSize.

После завершения цикла поле независимой переменной currt принимает конечное значение tEnd, делая последний шаг, и вызывается обработчик callBack.