Оптимальный режим управления двухсекторной моделью экономики

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Содержание

Введение

1. Постановка задачи

2. Задача оптимального управления. Достаточные условия

2.1 Постановка задачи оптимального управления

2.2 Принцип максимума Понтрягина

2.3 Достаточное условие Эрроу

3. Численное решение задачи

3.1 Метод Эйлера

3.2 Методы Рунге-Кутта III, IV порядков

3.3 Метод Адамса-Башфорта

4. Результат решения задачи оптимального управления

Список используемой литературы

Приложение А(обязательное). Текст программы



Введение

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

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

Математическому моделированию двухсекторной экономической модели посвящено немного работ. В основном только рассматривается аналитическое решение. В качестве примера можно привести работу такого исследователя как Andrea Calogero.

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

Объект исследования - двухсекторная экономическая модель (в чистом виде).

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

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

Для достижения поставленной цели необходимо решить следующие задачи:

1) Применить принцип максимума Понтрягина для решения задачи оптимального режима управления двухсекторной модели в экономике;

2) Применить достаточное условие для проверки полученной экстремали на оптимальность в задачи оптимального режима управления двухсекторной модели в экономике;

3) Представить непрерывную задачу оптимального режима управления двухсекторной модели в экономике в дискретном виде: Эйлера, Рунге-Кутта III, Рунге-Кутта IV, Адамса-Башфорта;

4) Разработать программное средство для численного решения задачи оптимального режима управления двухсекторной модели в экономике.



1. Постановка задачи

Рассмотрим двухсекторную экономическую модель, где сектор 1--производит инвестиционные товары, основные фонды, или капитал, фондообразующего сектора, осуществляющего инвестиции в собственное развитие и в развитие потребительского сектора (фазовая координата, заданная в начале свободная в конце, где она должна быть выбрана из соображений оптимальности);

сектор 2-- производит потребительские товары, основные фонды потребительского сектора, производящего товары

потребления (фазовая координата, отсчитываемая от достигнутого предпланового уровня);

-- доля инвестиций, направляемых в потребительский сектор (управление);

Т -- заданная протяженность интервала планирования, отсчитываемого от нуля.

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

(1.1)

с начальными условиями

, (1.2)

Где

- коэффициент степени воздействия на рост производства;

- объем производства в начальный момент времени инвестиционных товаров;

- объем производства в начальный момент времени потребительских товаров;

- производство в 1-ом секторе;

- производство во 2-ом секторе.

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

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

В частном случае, рассмотрим задачу максимизации объема потребления в данный плановый период [0, Т], при , тогда получим систему (1.3):

(1.3)



2. Задача оптимального управления. Достаточные условия

2.1 Постановка задачи оптимального управления

Требуется найти максимум функционала:

(2.1)

при динамических ограничениях:

(2.2)

(2.3)

начальных условиях:

, (2.4)

ограничениях на управление:

(2.5)

2.2 Принцип максимума Понтрягина

Составим Гамильтониан для задачи (2.1)-(2.5) (рассмотрим регулярный случай, то есть при ):

(2.6)

Так как на управление наложено ограничение, применим принцип максимума Понтрягина по теореме:

Рассмотрим и и пусть u^* - оптимальным управлением задачи (2.1)-(2.5), х^*- соответствующая ему оптимальная траектория. Тогда существуют множители, где - константа, непрерывные координаты, и выполняются условия:

1) принцип максимума Понтрягина (РМР)

Для управление в этой точке принадлежит , то есть Гамильтониан найдется в точке

(2.7)

2) сопряженнее уравнение

(2.8)

3) условие трансверсальности

(2.9)

4)

Тогда РМР для задачи (2.1)-(2.5) имеет вид:

(2.10)



2.3 Краевая задача принципа максимума

(2.11)

(2.12)

Так как задача (2.1)-(2.5) со свободным правым концом, условия трансверсальности необходимы

(2.13)

(2.14)

Найдем максимальное положение управления в гамильтониане Н из (2.6) во всей области значений переменных x₁, x₂ , . Рассмотрим зависимость Н по , зависимость линейна, следовательно, максимум достигается на концах отрезка.

Проинтегрируем (2.12), и, используя условие трансверсальности (2.13) найдем константу, тогда . Кроме того (2.13) и (2.14) в уравнение (2.11) дают . Заметим, что . Из непрерывности множителей следует, что существует такое, что

(2.15)

(2.16)

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

Найдем и

Проинтегрируем по t и получим следующее:

(2.17)

Вычисление показывает, что предположение (2.15) справедливо для . Поэтому давайте предположим, что существует точка такая, что

(2.18)

В силу (2.16) на интервале мы получаем, что . Теперь (2.11) дает, с учетом непрерывности в точке :

(2.19)

Легко проверить, что поскольку и , из выпуклости функции и следует, что предположение (2.18) выполняется с ф'= 0.

В результате управление удалось построить как функцию времени:

(2.20)

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

Сначала решим первое дифференциальное уравнение, поскольку оно не содержит , разделим переменные, проинтегрируем и выделим :

(2.21)

(2.22)

Подставив исходное уравнение и используя начальные условия, получим, что

(2.23)

(2.24)

Решая второе дифференциальное уравнение, проинтегрируем и выделим , подставив вместо полученное значение:

(2.25)

(2.26)

Используя динамику и начальное условие на траектории, получим решения подозрительные на оптимальность:

(2.27)

(2.28)

(2.29)

(2.30)

2.3 Достаточное условие Эрроу

Так как экстремальные решения являются лишь подозрительными на максимум, то необходимо проверить достаточные условия.

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

Рассмотрим задачу на максимум (2.1)-(2.5) с и , - должно быть нормально оптимальным управлением, х^* - соответствующая траектория и множители. Рассмотрим развернутую функцию Гамильтониана Н₀:

(2.31)

и предположим, что каждомуфункции

(2.32)

вогнута. Кроме того, мы предполагаем, что функция U вдоль кривой равен u^*, то есть

(2.33)

- оптимальное управление

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

(2.34)

u=(2.35)

Заметим, что для каждого фиксированных функции является вогнутой, и что функции совпадает с : достаточное условие Эрроу выполняется.



3. Численное решение

Будем рассматривать следующую задачу оптимального управления:

(3.1)

(3.2)

(3.3)

, (3.4)

, (3.5)

(3.6)

Для численного интегрирования задачи Коши (3.1) при условии существования и единственности решения, известны методы: одношаговые Эйлера и Рунге-Кутта с модификациями, многошаговые типа Адамса и др.

Рассмотрим алгоритм решения задачи оптимального управления.

Алгоритм

Шаг 1: Подается начальный вектор управлений- u, начальные условия -, область поиска - , точность вычислений -, начальный и конечный момент времени - [t,T].

Шаг 2: Выбираем точку переключения

Шаг 3: При данном управлении решается задача Коши для системы дифференциальных уравнений, ищется траектория, вычисляется значение целевого функционала.

Шаг 4: Вызывается функция поиска приближенного решения.

Шаг 5: Проверка условий останова. Если условия не выполняются,переходим к шагу 2,иначе алгоритм прекращает свою работу.



3.1 Метод Эйлера

Исторически первым и наиболее простым способом численного решения задачи Коши для ОДУ первого порядка является метод Эйлера. В его основе лежит аппроксимация производной отношением конечных приращений зависимой y и независимой x переменных между узлами равномерной сетки:

(3.1)

где y_i+1 это искомое значение функции в точке x_i+1.

Если теперь преобразовать это уравнение, и учесть равномерность сетки интегрирования, то получится итерационная формула, по которой можно вычислить y_i+1 , если известно y_i в точке х_i:

(3.2)

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

Построим рекуррентные формулы по формуле (3.2):

динамические ограничения:

(3.3)

(3.4)

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



(3.5)

условия трансверсальности:

(3.6)

ограничение на управление:

(3.7)

Гамильтониан Н:

, при (3.8)

Функция Лагранжа:

(3.9)

Условие стационарности по фазовым переменным:

(3.10)

(3.11)

(3.12)

(3.13)

Условие стационарности по управлению:

3.14)

3.2 Методы Рунге-Кутты III,IV порядков

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

Воспользовавшись формулой Симпсона, можно получить еще более точную формулу для решения задачи Коши ОДУ первого порядка - широко используемого в вычислительной практике метода Рунге-Кутты.

В формуле Симпсона для приближенного вычисления определенного интеграла используются значения подынтегрального выражения в трех точках. В интеграле их всего две, поэтому введем дополнительную точку в середине отрезка [x_i+1,x_i]

тогда можно переписать так:

Полученное выражение является неявным, так как в правой части содержатся еще не определенные значения функции y_i+h/2 и y_i+1. Чтобы воспользоваться этой формулой, надо использовать некоторое приближение для вычисления этих значений

(3.15)



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

Алгоритм Рунге-Кутты третьего порядка - РК3 (погрешность порядка h3):

(3.15)

Алгоритм Рунге-Кутты четвертого порядка- РК4(погрешность порядка h4):

(3.16)

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

Построим рекуррентные формулы по формуле по алгоритму РК3:

(3.10)

(3.11)

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

(3.12)

условия трансверсальности:

(3.13)

ограничение на управление:

(3.14)

Гамильтониан Н:

, при (3.15)

Условия стационарности по управлению и фазовым переменным строится аналогично методу Эйлера.

Для РК4 рекуррентные формулы строятся аналогично РК3, используя формулы (3.16).



3.3 Метод Адамса-Башфорта

Рассмотренные ранее методы (Эйлера, Рунге-Кутты) используют значение функции на одном предшествующем шаге, поэтому они относятся к так называемым одношаговым методам. Точность вычислений можно увеличить, если использовать при нахождении решения в некотором узле x_iинформацию о значениях функции, полученных в нескольких (k) предыдущих узлах сетки интегрирования (x_i-1, x_i-2… x_i-k).

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

(3.17)

Значения квадратурных коэффициентов для k от 2 до 4 приведены в таблице.

k
2	3/2	-1/2
3	23/12	-16/12	5/12
4	55/24	-59/24	37/24	-9/24

Полученное таким образом семейство формул называется явной k-шаговой схемой Адамса (методы Адамса-Башфорта).

Так как для вычислений по k-шаговой формуле необходимо знание значения функции в k узлах. Поэтому приходится (k-1) решение в первых узлах x₁, x₂, …, x_k-1получать с помощью какого-либо одношагового метода, например метода Рунге-Кутты 4-го порядка.



4. Результаты решения задачи оптимального управления

Для решения задачи оптимального управления на тему оптимальный режим управления двухсекторной модели в экономике разработано программное средство в среде Visual Studio 2013 на языке С# (приложение А). В качестве входных параметров будем рассматривать:

- коэффициент степени воздействия на рост производства;

- объем производства в начальный момент времени инвестиционных товаров;

- объем производства в начальный момент времени потребительских товаров;

- протяженность интервала планирования, отсчитываемого от нуля.

На рисунке (4.1) представлено решение задачи при следующих параметрах, методом Эйлера.

Рисунок 5.1 - Решение задачи оптимального управления методом Эйлера при .

Для примера на рисунке (5.2) представлено решение задачи оптимального управления методом более высокого порядка - метод Рунге-Кутты IV порядка.

Рисунок 5.2 - Решение задачи оптимального управления методом Рунге-Кутты IV при .

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

На рисунке (5.3) представлено сопоставление методов. Из реализованных методов наибольшую точность имеет метод Адамса-Башфорта.



Рисунок 5.3 - Сопоставление полученных решений численными методами.



Заключение

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

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

Было доказано с помощью достаточного условия Эрроу, что полученная экстремаль это и есть оптимальное решение.

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

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



Список использованных источников

1) CalogeroА. NOTES ON OPTIMAL CONTROL THEORY[текст].2 книга. - Universita di Milano-Bicocca, 2014 - 143 с.

2) Токарев В. В. Методы оптимальных решений. - М.: Изд-во ФИЗМАТЛИТ, 2011 - 422 с.

3) Фесько О. В. Программный комплекс оптимизации динамических систем на множествах управлений[текст].- Изд-во Вестник БГУ, 2010 -68 с.

4) Кротов В.Ф. Основы теории оптимального управления[текст]. - М.:Высш. Шк.,1990-430с.

5) Галеев Э.М. Курс лекций по вариационному исчислению и оптимальному управлению [текст]. - М.: МГУ.,1996 - 160 c.



Приложение А(обязательное)

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

using System;

using System.Collections.Generic;

using System.ComponentModel;

using System.Data;

using System.Drawing;

using System.Linq;

using System.Text;

using System.Threading.Tasks;

using System.Windows.Forms;

namespace WindowsFormsApplication4

{

public partial class Form1 : Form

{

public Form1()

{

InitializeComponent();

}

double x1_shtr(double t, double x1,double alfa,double u)

{

return x1 * alfa*u;

}

double x2_shtr(double t, double x1, double alfa, double u)

{

return x1 *(1-u)*alfa;

}

private void таблицаToolStripMenuItem_Click(object sender, EventArgs e)

{panel1.Visible = false;

label3.Text = "Решение представленно в таблице";

dataGridView1.Visible = true;

chart1.Visible = false;

chart2.Visible = false;

label1.Visible = false;

label2.Visible = false;

label14.Visible = false;

dataGridView2.Visible = false;

}

private void графикиToolStripMenuItem_Click(object sender, EventArgs e)

{

panel1.Visible = false;

label3.Visible = true;

label3.Text = "Решение представленно на графиках";

dataGridView1.Visible = false;

chart1.Visible = true;

chart2.Visible = true;

label1.Visible = true;

label2.Visible = true;

label5.Visible = true;

label14.Visible = false;

dataGridView2.Visible = false; }

private void

public void Out(int q, double[] t, double[] x1, double[] x2, double[] u1,double J1)

{int k = 0; int q1; q1 = q;

while (q1 != 0)

double[] u1 = new double[q]; ;

double[] psi1 = new double[q + 1];

double[] psi2 = new double[q + 1];

double J, J1;t[0] = 0;

for (int i = 0; i <= q - 1; i++)

t[i + 1] = t[i] + h;

for (int i = 0; i <= q - 1; i++)

u[i] = 0.5;

x1[0] = 1;

for (int i = 0; i <= q - 1; i++)

x1[i + 1] = x1[i] + (1 - u[i]) * x1[i] * h;

x2[0] = 0;

for (int i = 0; i <= q - 1; i++)

x2[i + 1] = x2[i] + u[i] * x1[i] * h;

J = -x2[q];label:

psi1[q] = 0; psi2[q] = 1;

for (int i = q - 1; i >= 1; i--)

{

psi1[i] = psi1[i + 1] + psi1[i + 1] * (1 - u[i]) * h + psi2[i + 1] * u[i] * h;

psi2[i] = psi2[i + 1];}

for (int i = 0; i <= q - 1; i++)

{if (psi1[i + 1] >= psi2[i + 1])

u1[i] = 0;

else u1[i] = 1;}

x1_1[0] = 1;

for (int i = 0; i <= q - 1; i++)

x1_1[i + 1] = x1_1[i] + (1 - u1[i]) * x1_1[i] * h;

x2_1[0] = 0;

for (int i = 0; i <= q - 1; i++)

x2_1[i + 1] = x2_1[i] + u1[i] * x1_1[i] * h;

J1 = -x2_1[q];

if (J1 < J){

if (Math.Abs(J1 - J) <= epsJ)

Out(q, t, x1_1, x2_1, u1,J1);

else{J = J1;

for (int i = 0; i <= q; i++)

{x1[i] = x1_1[i];

x2[i] = x2_1[i];}

for (int i = 0; i <= q - 1; i++)

u[i] = u1[i];

goto label; }

Else Out(q, t, x1, x2, u,J);}

private void button2_Click(object sender, EventArgs e)

{ double T;

T = Convert.ToDouble(textBox2.Text);

int q;q = 1000;double epsJ;

epsJ=0.0001;double h;

h = T / q;

double x1, x2, k1, k2, k3, k4, m1, m2, m3, m4, U;

double[] t = new double[q + 1];

double[] x1 = new double[q + 1];

double[] x2 = new double[q + 1];

double[] u = new double[q]; ;

double[] x1_1 = new double[q + 1];

double[] x2_1 = new double[q + 1];

double[] u1 = new double[q]; ;

double[] psi1 = new double[q + 1];

double[] psi2 = new double[q + 1];

double J, J1;t[0] = 0;

psi1[q] = 0; psi2[q] = 1;

for (int i = q - 1; i >= 1; i--)

{

psi1[i] = psi1[i + 1] + psi1[i + 1] *

(1 - u[i]) * h + psi2[i + 1] * u[i] * h;

psi2[i] = psi2[i + 1];}

for (int i = 0; i <= q - 1; i++)

{if (psi1[i + 1] >= psi2[i + 1])

u1[i] = 0;

else u1[i] = 1;}

x1_1[0] = 1;

for (int i = 0; i <= q - 1; i++)

{

k1 = x1_shtr(t, x1, alfa, U);

k2 = x1_shtr(t + tau / 2, x1 + tau / 4 * k1, alfa, U);

k3 = x1_shtr(t + tau / 2, x1 + tau / 2 * k2, alfa, U);

k4 = x1_shtr(t + tau, x1 + tau * k1 - 2 * tau * k2 + 2 *

tau * k3, alfa, U);

x1_1[i+1] = x1_1[i] + tau / 6 * (k1 + 2 * k2 + 2 * k3 + k4);}

t = t + tau;

t x2_1[0] = 0;

for (int i = 0; i <= q - 1; i++)

{

k1 = x1_shtr(t, x1, alfa, U);

k2 = x1_shtr(t + tau / 2, x1 + tau / 4 * k1, alfa, U);

k3 = x1_shtr(t + tau / 2, x1 + tau / 2 * k2, alfa, U);

k4 = x1_shtr(t + tau, x1 + tau * k1 - 2 * tau * k2 + 2 *

tau * k3, alfa, U);

x2_1[i+1] = x2_1[i] + tau / 6 * (k1 + 2 * k2 + 2 * k3 + k4);

}

J1 = -x2_1[q];

if (J1 < J){

if (Math.Abs(J1 - J) <= epsJ)

Out(q, t, x1_1, x2_1, u1,J1);

else{J = J1;

for (int i = 0; i <= q; i++)

{x1[i] = x1_1[i];

x2[i] = x2_1[i];}

for (int i = 0; i <= q - 1; i++)

u[i] = u1[i];

goto label; }

Else Out(q, t, x1, x2, u,J);

//ТОЧНОЕ РЕШЕНИЕ 1

t = 0; k = 0;

while (t < T){if ((t >= 0) && (t <= (T - 2 / alfa)))

{x1_shtr_tochnoe = alfa1 * Math.Exp(alfa * t);

dataGridView1.Rows[k].Cells[2].Value =

x1_shtr_tochnoe;

t = t + tau;

k++; chart1.Series[1].Points.AddXY(t,

x1_shtr_tochnoe);}

else{

x1_shtr_tochnoe = alfa1 * Math.Exp(alfa * T - 2);

dataGridView1.Rows[k].Cells[2].Value =

x1_shtr_tochnoe;

t = t + tau;k++;

chart1.Series[1].Points.AddXY(t,

x1_shtr_tochnoe);}}

//ТОЧНОЕ РЕШЕНИЕ 2t = 0; {q1 = q1 / 10; k++;}

chart1.Series.Clear();

chart2.Series.Clear()

chart1.Series.Add("x1(t)");

chart1.Series.Add("x2(t)");

chart2.Series.Add("u(t)"); chart1.Series["x1(t)"].ChartType = SeriesChartType.Line; chart1.Series["x2(t)"].ChartType = SeriesChartType.Line;

chart2.Series["u(t)"].ChartType = SeriesChartType.Line;

chart1.Series["x1(t)"].BorderWidth = 2;

chart1.Series["x2(t)"].BorderWidth = 2;

chart2.Series["u(t)"].BorderWidth = 2;

for (int i = 0; i <= q; i++){

chart1.Series["x1(t)"].Points.AddXY(t[i], x1[i]); chart1.Series["x2(t)"].Points.AddXY(t[i], x2[i]); }

for (int i = 0; i <= q - 1; i++)

chart2.Series["u(t)"].Points.AddXY(t[i], u1[i]);

textBox1.Text = Convert.ToString(J1);

string[] row = new string[5];

row[0] = "q";

row[1] = "t"; row[2] = "x(t)";

row[3] = "y(t)";

row[4] = "u(t)";

dataGridView1.ColumnCount = 5;

dataGridView1.RowCount = 1;

dataGridView1.Rows.Add(row);

for (int i = 0; i <= q; i++){

row[0] = i.ToString();

row[1] = Math.Round(t[i], k).ToString();

row[2] = x1[i].ToString();

row[3] = x2[i].ToString();

if (i != q)

{row[4] = u1[i].ToString();}

else

{row[4] = "";}

dataGridView1.Rows.Add(row);}}

double x1_shtr(double t, double x1,double alfa,double u)

{ return x1 * alfa*u}

double x2_shtr(double t, double x1, double alfa, double u)

{return x1 *(1-u)*alfa;}

private void button1_Click(object sender, EventArgs e)

{ double T;

T = Convert.ToDouble(textBox2.Text);

int q;q = 1000;double epsJ;

epsJ=0.0001;double h;

h = T / q;

double[] t = new double[q + 1];

double[] x1 = new double[q + 1];

double[] x2 = new double[q + 1];

double[] u = new double[q]; ;

double[] x1_1 = new double[q + 1];

double[] x2_1 = new double[q + 1];

for (int i = 0; i <= q - 1; i++)

t[i + 1] = t[i] + h;

for (int i = 0; i <= q - 1; i++)

u[i] = 0.5;

x1[0] = 1;

for (int i = 0; i <= q - 1; i++)

{

k1 = x1_shtr(t, x1[i], alfa, U[i]);

k2 = x1_shtr(t + tau / 3, x1[i] + tau / 3 * k1, alfa, U[i]);

k3 = x1_shtr(t + 2*tau / 3, x1[i] + 2*tau / 3 * k2, alfa, U[i]);

x1[i+1] = x1[i] + tau / 4 * (k1 + 3 * k3);

t = t + tau;

}

x2[0] = 0;

for (int i = 0; i <= q - 1; i++)

{

k1 = x2_shtr(t, x1[i], alfa, U[i]);

k2 = x2_shtr(t + tau / 3, x1[i] + tau / 3 * k1, alfa, U[i]);

k3 = x2_shtr(t + 2*tau / 3, x1[i] + 2*tau / 3 * k2, alfa, U[i]);

x2[i+1] = x2[i] + tau / 4 * (k1 + 3 * k3);

}

J = -x2[q];label:

psi1[q] = 0; psi2[q] = 1;

for (int i = q - 1; i >= 1; i--)

{

psi1[i] = psi1[i + 1] + psi1[i + 1] * (1 - u[i]) * h + psi2[i + 1] * u[i] * h;

psi2[i] = psi2[i + 1];}

for (int i = 0; i <= q - 1; i++)

{if (psi1[i + 1] >= psi2[i + 1])

u1[i] = 0;

else u1[i] = 1;}

x1_1[0] = 1;

for (int i = 0; i <= q - 1; i++)

{

k1 = x1_shtr(t, x1[i]_1, alfa, U[i]);

k2 = x1_shtr(t + tau / 3, x1[i]_1 + tau / 3 * k1, alfa, U[i]);

k3 = x1_shtr(t + 2*tau / 3, x1[i]_1 + 2*tau / 3 * k2, alfa, U[i]);

x1_1[i+1] = x1_1[i] + tau / 4 * (k1 + 3 * k3);

t x2_1[0] = 0;

for (int i = 0; i <= q - 1; i++)

{

k1 = x2_shtr(t, x1[i], alfa, U[i]);

k2 = x2_shtr(t + tau / 3, x1[i] + tau / 3 * k1, alfa, U[i]);

k3 = x2_shtr(t + 2*tau / 3, x1[i] + 2*tau / 3 * k2, alfa, U[i]);

x2_1[i+1] = x2_1[i] + tau / 4 * (k1 + 3 * k3);

}

J1 = -x2_1[q];

if (J1 < J){

if (Math.Abs(J1 - J) <= epsJ)

Out(q, t, x1_1, x2_1, u1,J1);

else{J = J1;

for (int i = 0; i <= q; i++)

{x1[i] = x1_1[i];

x2[i] = x2_1[i];}

for (int i = 0; i <= q - 1; i++)

u[i] = u1[i];

label; }

x2[i+1] = x2[i] + tau / 6 * (k1 + 2 * k2 + 2 * k3 + k4);

t = t + tau;

}

J = -x2[q];label:

k = 0;

//погрешность между точным и

численным решением х1*

double max = Math.Abs(Convert.

ToDouble(dataGridView1.Rows[0].Cells[2].Value)

- Convert.ToDouble(dataGridView1.Rows[0].Cells[0].

Value));

for (int i = 1; i < n; i++){if (max < (Math.Abs(Convert

.ToDouble(dataGridView1.Rows[i]}

Размещено на Allbest.ru