Добавить материал и получить бесплатное свидетельство о публикации в СМИ
Эл. №ФС77-60625 от 20.01.2015

Опубликуйте свой материал в официальном Печатном сборнике методических разработок проекта «Инфоурок»

(с присвоением ISBN)

Выберите любой материал на Вашем учительском сайте или загрузите новый

Оформите заявку на публикацию в сборник(займет не более 3 минут)

+

Получите свой экземпляр сборника и свидетельство о публикации в нем

Инфоурок / Информатика / Конспекты / Лабораторный практикум по моделированию
ВНИМАНИЮ ВСЕХ УЧИТЕЛЕЙ: согласно Федеральному закону № 313-ФЗ все педагоги должны пройти обучение навыкам оказания первой помощи.

Дистанционный курс "Оказание первой помощи детям и взрослым" от проекта "Инфоурок" даёт Вам возможность привести свои знания в соответствие с требованиями закона и получить удостоверение о повышении квалификации установленного образца (180 часов). Начало обучения новой группы: 24 мая.

Подать заявку на курс
  • Информатика

Лабораторный практикум по моделированию

библиотека
материалов


hello_html_m53d4ecad.gifhello_html_m53d4ecad.gifМИНОБРНАУКИ РОССИИ

Федеральное государственное бюджетное образовательное

учреждение высшего профессионального образования

«Чувашский государственный университет имени И. Н. Ульянова»









В. К. Никишев





Математическое

моделирование



Лабораторный практикум






Чебоксары

2013


УДК 004.92(076.5)

ББК 3973.2-044.4я73

Н62

.



Никишев В. К.

Н62 Математическое моделирование: лабораторный практикум.

Чебоксары: Изд-во Чуваш. Ун-та, 2013. – 151 с.



Представлены примеры выполнения лабораторных работ по моделированию на языках программирования Delphi, VS 2010 ( VC#,VC++, VB.NET) и в среде моделирования SciLab (MatLab). Тематика лабораторных работ соответствует рабочей программе по математическому моделированию, которая составлена в соответствии со стандартом образования. Каждая работа содержит условие задачи, алгоритм в виде блок-схемы, программы на языках программирования VC#, VC++, DELPHI и результаты исследований. Задания для выполнения лабораторных работ представлены в конце лабораторного практикума.

Для бакалавров II-III, магистров и аспирантов технических факультетов, изучающих математическое моделирование с использованием современных языков.



Ответственный редактор канд. техн. наук, профессор В. К. Никишев


Утверждено Учебно-методическим советом университета

УДК 004/92(076/5)


ISBN 978-5-7677-1739-2 © Издательство Чувашского

Университета, 2013 © Никишев В. К., 2013






Предисловие

В настоящее время большое внимание уделяется вопросам моделирования различных систем с использованием современных языков программирования (Visual Basic, Delphi, VC#, VC++, VB.NET) и информационных программ, например, Excel, MathCad, MatLab, Maple, SciLab.

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












Лабораторная работа 1. Методы исследования объектов, динамика которых описывается дифференциальными уравнениями с использованием языков программирования

Delphi, VC++.NET, VC#. NET

Цель занятия:

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

2. Научиться разрабатывать алгоритм и программу с использованием языков программирования Delphi, VC++.NET, VC#.NET

3. Практически усвоить численные методы Эйлера и Рунге-Кутта для решения дифференциальных уравнений 1-го порядка.

Задачи занятия:

1. Разработка алгоритма в виде блок-схемы.

2. Построение графиков кривых y=f(x), ý = f(x) при параметрах a-const и var.

3. Анализ результатов исследований.

Модели объектов исследования

aӳ+bý + cy=f

Программа исследования

1. a,c,f - const, t – var ( t0 – tk, h= 0.1, 0.01)

2. a,f - const, t -var( t0 – tk, h= 0.1, 0.01), b- var

Пример 1

Условие задачи: cоставить алгоритм и проект моделирования объекта, динамика которого описывается дифференциальным уравнением 1-го порядка методом Эйлера

5 ý + 3 y = 4/








Блок-схема алгоритма



hello_html_10fa2fed.png




































Отчет по лабораторной работе


hello_html_m66970af2.pnghello_html_197f5ab1.png


Рис.1. Форма 2. Автор Рис. 2. Форма 4. Условие


hello_html_5682e6ae.png








Рис. 3. Форма 1. Титульный лист


Форма для исследования объекта

hello_html_m4ddae869.png













Рис. 4. Форма 3. Результаты исследования


Листинг программы ( язык программирования Delphi)

procedure TForm3.Button2Click(Sender: TObject); ( рис. 4).

var: x:real; i:integer;

function F(x,y:real):real;

begin

f:=5*x+(3*y)-4;

end;

begin

a:=StrToFloat(Edit1.text);

b:=StrToFloat(Edit2.text);

h:=StrToFloat(Edit3.text);

y1:=strtofloat(edit4.text); //шаг

memo1.lines.add(floattostr(x));

memo2.lines.add(floattostr(y));

x:=a;

Chart1.Series[0].clear; Chart1.Series[1].clear;

repeat

begin

y:=y1+h*f(x,y);

Chart1.Series[0].Add(y,FloatToStr(x),clRed);

Chart1.Series[1].Add(f(x,y),FloatToStr(x),clblue);

x:=x+h; y1:=y;

memo1.lines.add(floattostrf(x,fffixed,7,3));

memo2.lines.add(floattostrf(y,fffixed,7,3));

end;

until x>b;

end;


Листинг программы ( язык программирования VC++)


Исходное дифференциальное уравнение


ý - y - 2*sin(x) = 0;


hello_html_647c1997.png















Рис 5. Форма для условия задачи.


private: System::Void button1_Click(System::Object^ sender, System::EventArgs^ e) {

this->Hide();

}

private: System::Void button2_Click(System::Object^ sender, System::EventArgs^ e) {

double f1,x, y1;

double h = Convert::ToDouble(textBox1->Text);

double b = Convert::ToDouble(textBox2->Text);

double y = Convert::ToDouble(textBox3->Text);

x = 0;

y1 = y;

chart1->Series[0]->Points->Clear();

chart1->Series[1]->Points->Clear();

do

{

f1 = y - 2*sin(x);

x = x+h;

y = y+f1*h;

chart1->Series[0]->Points->AddXY(x, y);

chart1->Series[1]->Points->AddXY(x,f1);

}

hello_html_3edb494b.png while (x<=b);

}};}};}










Рис 6. Форма для автора.


hello_html_17802094.png










Рис 7. Форма для отображения результатов моделировния.

Язык программирования C#

private void button1_Click(object sender, EventArgs e)

{

double h = 0.01;

double y = 0;

double b = 2;

double x = 0;

double y1 = 4;

do

{

double f = y + 2*Math.Sin(x);

x = x + h;

y = y1 + f * h;

chart1.Series[0].Points.AddXY(x, y);

chart1.Series[1].Points.AddXY(x, f);

y1 = y;

}

while (x <= b);

}

// foreach (DataPoint p in chart1.Series[0].Points)

// {

// Вывожу Х в лог

// textBox1.AppendText("X=" + p.XValue.ToString());

// textBox1.AppendText(Environment.NewLine);

// listBox1.Items.Add("X=" + p.XValue.ToString());

// }

// Y является массивом, поэтому пробегаю по массиву

// foreach (DataPoint yp in chart1.Series[0].Points)

// {

//Вывожу Y

// textBox1.AppendText("Y=" + yp.ToString());

// textBox1.AppendText(Environment.NewLine);

// listBox1.AppendText("Y=" + yp.ToString());

// listBox1.AppendText(Environment.NewLine);

// listBox2.Items.Add("Y=" + yp.ToString());

// }

// }

private void chart1_Click(object sender, EventArgs e)

{

}

private void button2_Click_1(object sender, EventArgs e)

{

Close();

} }}







Пример 2.

Условия задачи: Исследовать объект, динамика которого представляется дифференциальным уравнением


ý - y - a*sin(x) = 0;


Цель занятия: получить практические навыки исследования систем (объектов), динамика которых описывается дифференциальными уравнениями 1-го порядка при изменении параметра «а»

Задачи занятия:

1. Разработка алгоритма в виде блок- схемы.

2. Построение графиков кривых y=f(x), ý = f(x) =y(x) при параметрах a-var


hello_html_5d754266.png

hello_html_m71f2d463.png












Рис 8. Титульный лист Рис. 9. Условия задачи


hello_html_789953ac.png



Рис 10. Автор программы




hello_html_m377ae571.png















Рис. 11. Результаты исследования.

Листинг программы

private: System::Void button1_Click(System::Object^ sender, System::EventArgs^ e) {

this->Hide();

}

private: System::Void button2_Click(System::Object^ sender, System::EventArgs^ e) {

double f1,f2,x, y1;

double h = Convert::ToDouble(textBox1->Text);

double b = Convert::ToDouble(textBox2->Text);

double y2 = Convert::ToDouble(textBox3->Text);

double a = Convert::ToDouble(textBox4->Text);

double a1 = Convert::ToDouble(textBox5->Text);

x = 0;

// double y1 = y2;

int a3=a;

//for int i=0;i<=1;i++)

//{

chart1->Series[0]->Points->Clear();

chart1->Series[1]->Points->Clear();

chart1->Series[2]->Points->Clear();

chart2->Series[0]->Points->Clear();

chart2->Series[1]->Points->Clear();

chart2->Series[2]->Points->Clear();

//}

do

{

// for (int i=1;i<=5;i++)

// {

a3=a;

f1 = y2 - a3*sin(x);

// int i=0;

double y = y2+f1*h;

//int i=0;

chart1->Series[0]->Points->AddXY(x, y);

chart2->Series[0]->Points->AddXY(x,f1);

a3=a3+a1;

f1 = y2 - a3*sin(x);

y = y2+f1*h;

chart1->Series[1]->Points->AddXY(x,y);

chart2->Series[1]->Points->AddXY(x,f1);

a3=a3+a1;

f1 = y2 - a3*sin(x);

y = y2+f1*h;

chart1->Series[2]->Points->AddXY(x,y);

chart2->Series[2]->Points->AddXY(x,f1);

//y=y2;

x = x+h;

}

while (x<=b);

x=0;

}

private: System::Void chart2_Click(System::Object^ sender, System::EventArgs^ e) {

}

private: System::Void textBox1_TextChanged(System::Object^ sender, System::EventArgs^ e) {

}

};

}

Лабораторная работа 2

Методы исследования объектов, динамика которых описывается дифференциальными уравнениями 2-го поряда.

Цель занятия: Получить практические навыки исследования систем (объектов), динамика которых описывается дифференциальными уравнениями 2-го порядка.

Задачи занятия:

1. Разработка алгоритма в виде блок - схемы

2. Построение графиков кривых y=f(x), ý = f(x) =y(x) при параметрах a-const,

3. Анализ результатов исследований.


Пример исследования на языке программирования VC++


hello_html_m3e94ca83.png











Рис. 12. Титульный лист


hello_html_m574ca861.pnghello_html_6fac151e.png










Рис 13. Автор программы Рис 14. Условие задачи

hello_html_m887d5c3.png










Лист 15. Результаты исследования



Листинг программы

private: System::Void button1_Click(System::Object^ sender, System::EventArgs^ e) {

this->Hide();

}

private: System::Void button2_Click(System::Object^ sender, System::EventArgs^ e) {

// int i = 1;

// int n = 100;

double f1,x,y11,y22;

double h = Convert::ToDouble(textBox1->Text);

double b = Convert::ToDouble(textBox2->Text);

double y1 = Convert::ToDouble(textBox3->Text);

double y2 = Convert::ToDouble(textBox4->Text);

x = 0;

chart1->Series[0]->Points->Clear();

chart1->Series[1]->Points->Clear();

do

{

// f1 = y - 2*sin(x);

x = x+h;

// y = y+f1*h;

y11=y1+h*y2;

y22=y2+h*(5-3*y2-4*y1);

//y = y1+h*(f1+(y - 2*x/y))/2;

chart1->Series[0]->Points->AddXY(x, y11);

chart1->Series[1]->Points->AddXY(x,y22);

y1=y11;

y2=y22;

}

while (x<=b);

}

private: System::Void Form5_Load(System::Object^ sender, System::EventArgs^ e) {

}};}


Пример исследования на языке Delphi

Условие задачи: cоставить алгоритм и проект моделирования объекта, динамика которого описывается дифференциальным уравнением 2-го порядка методом Рунге-Кутта

y’’+a*x*y’-y=0,4.

hello_html_m61431df2.gif









hello_html_481adf70.gif










Рис. 15. Титульный лист и условие задачи


hello_html_d51a769.gif











hello_html_m36fa6c09.gif



Рис. 16. Задачи моделирования






















Рис. 17. Результаты моделирования




Блок-схема алгоритма




hello_html_5da2e18e.png

Программа на Delphi

unit Unit2;

interface

var

Form2: TForm2;

implementation

uses Unit1,Unit4;

procedure TForm2.SpeedButton1Click(Sender: TObject);

begin

form4.Show;form2.Hide;

end;

procedure TForm2.SpeedButton2Click(Sender: TObject);

var

i,p :integer; // переменные, использующиеся в циклах

a : integer; // а - параметр а в нашем уравнении

h,t:real; // h - точность, t - шаг

n :integer; // число шагов

// x0,y0 :real; // начальные значения x и y

tn, tk : integer; // границы

x,y,kx1,kx2,kx3,kx4,ky1,ky2,ky3,ky4,dx,dy :array[1..10000] of real;

// x,y - значения x и y

// kx1-4,ky1-4 - переменные, принимающие участие в формуле Рунге-Кутта

//dx, dy - отрезки x и y решение уравнения по Рунге-Кутта

begin

a:=strtoint(edit1.Text);tn:=strtoint(edit2.Text);

tk:=strtoint(edit3.Text);h:=strtofloat(edit4.Text);

x[1]:=strtofloat(edit5.Text);y[1]:=strtofloat(edit6.Text);

n:=trunc((tk-tn)/h);

t:=tn;

chart1.SeriesList[0].Clear;chart2.SeriesList[0].Clear;

stringgrid1.RowCount:=n;

for i:= 1 to (n) do

begin

kx1[i]:=h*(-a*y[i]-6*x[i]);

kx2[i]:=h*(-a*y[i]-6*(x[i]+kx1[i]/2));

kx3[i]:=h*(-a*y[i]-6*(x[i]+kx2[i]/2));

kx4[i]:=h*(-a*y[i]-6*(x[i]+kx3[i]));

dx[i]:=1/6*(kx1[i]+2*kx2[i]+2*kx3[i]+kx4[i]);

x[i+1]:=x[i]+dx[i];

ky1[i]:=h*x[i];

ky2[i]:=h*(x[i]+ky1[i]/2);

ky3[i]:=h*(x[i]+ky2[i]/2);

ky4[i]:=h*(x[i]+ky3[i]);

dy[i]:=1/6*(ky1[i]+2*ky2[i]+2*ky3[i]+ky4[i]);

y[i+1]:=y[i]+dy[i];

t:=t+h;

stringgrid1.Cells[0,i]:=floattostr(t);

stringgrid1.Cells[1,i]:=floattostr(y[i]);

stringgrid1.Cells[2,i]:=floattostr(x[i]);

chart1.SeriesList[0].AddXY(t,y[i],'',clblue);

chart2.SeriesList[0].AddXY(t,x[i],'',clblue);

end;end;

procedure TForm2.FormCreate(Sender: TObject);

begin

stringgrid1.Cells[0,0]:='x';

stringgrid1.Cells[1,0]:='f(x)';

stringgrid1.Cells[2,0]:='f`(x)';end;end.

unit Unit5;

interface

var

Form5: TForm5;

implementation

uses Unit2, Unit4;

procedure TForm5.SpeedButton2Click(Sender: TObject);

var

i,p :integer; // переменные, использующиеся в циклах

a,a1,da : real; // а - параметр а в нашем уравнении

h,t:real; // h - точность, t - шаг

n,j :integer; // число шагов

// x0,y0 :real; // начальные значения x и y

tn, tk : integer; // границы

x,y,kx1,kx2,kx3,kx4,ky1,ky2,ky3,ky4,dx,dy :array[1..10000] of real;

// x,y - значения x и y на соотвествующем шаге

// kx1-4,ky1-4 - переменные опринимающие участие в формуле Рунге-Кутта

//dx, dy - отрезки x и y решение уравнения

begin

for j:= 0 to 5 do begin

chart1.SeriesList[j].Clear;chart2.SeriesList[j].Clear;

end;

a:=strtofloat(edit1.Text);a1:=strtofloat(edit7.Text);

da:=strtofloat(edit8.Text);

j:=0;a:=a-da;

repeat

tn:=strtoint(edit2.Text);tk:=strtoint(edit3.Text);

h:=strtofloat(edit4.Text);x[1]:=strtofloat(edit5.Text);

y[1]:=strtofloat(edit6.Text);n:=trunc((tk-tn)/h);

t:=tn;chart1.SeriesList[0].Clear;chart2.SeriesList[0].Clear;

for i:= 1 to (n) do

begin

kx1[i]:=h*(-a*y[i]-6*x[i]);

kx2[i]:=h*(-a*y[i]-6*(x[i]+kx1[i]/2));

kx3[i]:=h*(-a*y[i]-6*(x[i]+kx2[i]/2));

kx4[i]:=h*(-a*y[i]-6*(x[i]+kx3[i]));

dx[i]:=1/6*(kx1[i]+2*kx2[i]+2*kx3[i]+kx4[i]);

x[i+1]:=x[i]+dx[i]; ky1[i]:=h*x[i];

ky2[i]:=h*(x[i]+ky1[i]/2); ky3[i]:=h*(x[i]+ky2[i]/2);

ky4[i]:=h*(x[i]+ky3[i]);

dy[i]:=1/6*(ky1[i]+2*ky2[i]+2*ky3[i]+ky4[i]);

y[i+1]:=y[i]+dy[i]; t:=t+h;

chart1.SeriesList[j].AddXY(t,y[i],'',clblue);

chart2.SeriesList[j].AddXY(t,x[i],'',clblue);

end;

a:=a+da; j:=j+1;

until a > a1

end;

procedure TForm5.SpeedButton1Click(Sender: TObject);

begin

form4.show;form5.Hide;

end;end.


unit Unit6;

interface

var i,p :integer; // переменные, использующиеся в циклах

h,t:real; // h - точность, t - шаг

n :integer; // число шагов

// x0,y0 :real; // начальные значения x и y

tn, tk : integer; // границы

x,y,kx1,kx2,kx3,kx4,ky1,ky2,ky3,ky4,dx,dy: array[1..10000] of real;

// x,y - значения x и y на соответствующем шаге

// kx1-4,ky1-4 - переменные принимающие участие в формуле Рунге-Кутта

//dx, dy - отрезки x и y, из которых складывается решение уравнения по Рунге-Кутту

a: real; // a - parametr

begin

tn:=strtoint(edit2.Text);tk:=strtoint(edit3.Text);

h:=strtofloat(edit4.Text);x[1]:=strtofloat(edit5.Text);

y[1]:=strtofloat(edit6.Text);n:=trunc((tk-tn)/h);

t:=tn;

chart1.SeriesList[0].Clear;chart2.SeriesList[0].Clear;

a:=4;

for i:= 1 to (n) do

begin

kx1[i]:=h*(-a*y[i]-6*x[i]);

kx2[i]:=h*(-a*y[i]-6*(x[i]+kx1[i]/2));

kx3[i]:=h*(-a*y[i]-6*(x[i]+kx2[i]/2));

kx4[i]:=h*(-a*y[i]-6*(x[i]+kx3[i]));

dx[i]:=1/6*(kx1[i]+2*kx2[i]+2*kx3[i]+kx4[i]);

x[i+1]:=x[i]+dx[i]; ky1[i]:=h*x[i];

ky2[i]:=h*(x[i]+ky1[i]/2);

ky3[i]:=h*(x[i]+ky2[i]/2);

ky4[i]:=h*(x[i]+ky3[i]);

dy[i]:=1/6*(ky1[i]+2*ky2[i]+2*ky3[i]+ky4[i]);

y[i+1]:=y[i]+dy[i];

t:=t+h;

chart1.SeriesList[0].AddXY(t,y[i],'',clblue);

chart2.SeriesList[0].AddXY(t,x[i],'',clblue); end;end;


Лабораторная работа 3

Методы исследования объектов, динамика которых описывается дифференциальными уравнениями с использованием программы для моделирования

SciLab (MatLab)

Цель занятия:

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

2. Научиться разрабатывать алгоритм и программу с среде моделирования SciLab ( MatLab).

3. Практически усвоить численные методы Эйлера для решения дифференциальных уравнений 1-го и 2-го порядков

Задачи занятия:

1. Разработка алгоритма в виде блок-схемы

2. Построение графиков кривых y(x), dy/dx при параметрах a-const и var.

3. Анализ результатов исследований.


Модели объектов исследования


aӳ+bý + cy=f

Программа исследования

1. a,c,f - const, t – var ( t0 – tk, h= 0.1, 0.01)

2. a,f - const, t -var( t0 – tk, h= 0.1, 0.01), b- var


Пример 1

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

ý + 4 ysin(t)-5.


программа

function yd=f(t, y),yd=5*t-4*y*sin(t),endfunction;

y0=5;t0=0;t=0:0.01:3;

y=ode(y0,t0,t,f);

plot(t,y)

hello_html_m53bf2595.png

Рис. 20. Результаты исследования


Пример 2

Условие задачи: составить алгоритм и проект моделирования объекта, динамика которого описывается дифференциальным уравнением 1-го порядка методом Эйлера

ý + а* ysin(t)-5 , при изменении параметра на 4 значения.

N=4;

//disp(’Vvod N’);

//Цикл для ввода элементов в массиве y.

a=1;a1=2;

for i=1:N

function yd=f(t, y), yd=5*t-a*y*sin(t), endfunction;

y0=5;t0=0;t=0:0.01:3;

y=ode(y0,t0,t,f);

plot(t,y)

a=a+a1;

end

disp(y);




hello_html_1e20dc6b.png

Рис. 21. Результаты исследования



Пример 3

Условие задачи: составить алгоритм и проект моделирования объекта, динамика которого описывается дифференциальным уравнением 2-го порядка методом Эйлера

ӱ + 0.2*ý+ 4*y=2.


function dy=syst(t, y)dy=zeros(2,1);

dy(1)=y(2);

dy(2)=2-0.2*y(2)-4*y(1)

endfunction

y0=[0;0];t0=0;t=0:0.1:5;

y=ode(y0,t0,t,syst);

plot(t,y)

//end

//disp(y);









hello_html_m20de17d1.png












Рис. 21. Результаты исследования

Пример 4

Условия задачи. Исследовать объект, динамика которого описывается дифференциальным уравнением

ӱ + 0.2*ý+ 4*y=2

методом визуального моделирования

Исходное уравнение преобразуем к виду


Ӱ=2-0.2*ý-4*y

hello_html_m5efc0209.png














Рис. 22. Визуальная схема моделирования

hello_html_1a5ccb82.png

















Рис. 23. Результаты моделирования



hello_html_41096165.png














Рис. 24. Визуальная схема моделирования







hello_html_4d42c450.png



















Рис. 25. Результаты моделирования



Пример 5

Условие задачи: составить алгоритм и проект моделирования объекта, динамика которого описывается дифференциальным уравнением 2-го порядка методом Эйлера

ӱ + 0.2*ý+ 4*y=2


N=4;

//disp(’Vvod N’);

//Цикл для ввода элементов в массиве y.

//a=1;a1=5;

//for i=1:N

function dy=syst(t, y)

dy=zeros(2,1);

dy(1)=y(2);

dy(2)=2-0.2*y(2)-4*y(1)

endfunction

y0=[0;0];t0=0;t=0:0.1:5;

y=ode(y0,t0,t,syst);

plot(t,y)

//end

hello_html_m20de17d1.png



















Рис. 25. Результаты моделирования


Пример

Исследование в среде MatLab ( Simulink).

Условие задачи: составить алгоритм и проект:

моделирования объекта, динамика которого описывается

дифференциальным уравнением 2 порядка

y’’+axy’-y=0,4.




hello_html_2ac88e07.gif

2.2



























Рис 26. Визуальная схема и результаты исследования


Лабораторная работа 4. Исследования физических, биологических и других систем

Пример исследования физических систем. Свободное падение тел.

Условие задачи.

Цель занятия: изучить движение тел с высоты h при различных начальных условий: коэффициентов k1,k2, r

( радиуса шарика), плотности среды и формы тела.

Задачи занятия:

1. Получить зависимости h=f(t), v=f(t) при начальных значениях r,v0,h0,k1,k2.

2. Получить семейство графиков h=f(t,r), v=f(t,r).

3. Получить зависимость h=f(r), v=f(r) при t – const.

Модель объекта

hello_html_m53d4ecad.gifhello_html_m53d4ecad.gifhello_html_m76b23b0d.gif


Параметры k1 и k2 зависят от свойств среды и геометрии тела. Например, для шарика k1 =6μπr (формула Стокса), где μ– динамическая вязкость среды, r радиус шарика. Так, для воздуха μ = 0,0182 Н*с*м-2 при t=20o C и давлении 1 атм. Соответственно для воды μ = 1,002 Н*с*м-2. Для малых скоростей составляющей k1 можно пренебречь и = k2*v2, где k2 = 0,5с ; ρ плотность среды ; с коэффициент лобового сопротивления, который зависит от формы тела. Так, для диска в форме прямоугольника с=1,11, для сферы с=1,33 и с=0,55 в зависимости от направления расположения сферы, для шара с=0,4, а для каплевидных тел с=0,045 с=0,1 в зависимости от направления тела.






.





hello_html_14eded02.png















Пример. Исследовать падение шарика радиуса R с высоты h

Задачи исследования:

1. Получить зависимости h=f(t), v=f(t) при начальных значениях r,v0,h0,k1,k2.

2. Получить семейство графиков h=f(t,r), v=f(t,r).

3. Получить зависимость h=f(r), v=f(r) при t - const

Исходные данные:

  1. Модель объекта падения, R - радиус шарика, h - высота падения, v0 - начальная скорость падения, k1 k2 - коэффициенты внешней среды.


hello_html_6fdde3f.gif







  1. - площадь сечения тела, поперечного по отношению к потоку

  2. hello_html_4983766e.gif плотность шарика;

  3. hello_html_13faa.gif– плотность воздуха.

  4. c = 0.4 – коэффициент лобового сопротивления шарика;




hello_html_3c10f23b.gif

































var

Form4: TForm4;

km,r,m,vo,ho,s,p,v,a,dv,da,c:double ;

h1,h2,cn,ck,dc,t,dt,y1,y2,dy1,dy2:double;

ind,k:integer;

arrt,arrv,arra,arrh:array[1..10000] of double;

implementation

uses Unit1, Unit5;

//описание функции

function ff(m,km,r,c,s,p,v:double):double;

begin

ff:=(m*9.8-6*km*Pi*r*v-0.5*c*s*p*v*v)/m;

end;

//Кнопка возврата к титульному листу

procedure TForm4.N3Click(Sender: TObject);

begin

form5.show;

for ind:=0 to 7 do

begin

form5.Chart1.SeriesList[ind].Clear;

form5.Chart2.SeriesList[ind].Clear;

end;

//Входные параметры

r:=strtofloat(edit1.Text);//радиус падующего тела

km:=strtofloat(edit2.text);//коэффициент вязкости

s:=strtofloat(edit3.Text);//площадь падующего тела

p:=strtofloat(edit4.Text);//плотность среды

ho:=strtofloat(edit5.Text);//начальная высота

vo:=strtofloat(edit6.Text);//начальная скорость

m:=strtofloat(edit7.Text);//масса падающего тела

c:=0.1;//коэффициент лобового сопротивления среды

k:=0;//номер графика

// 1

while c<=1.5 do

begin

ind:=0;

t:=0;dt:=0.1;

y1:=vo;

y2:=ff(m,km,r,c,s,p,vo);

while t<=10 do

begin

form5.Chart1.SeriesList[k].Add(y1,floattostr(t));

dy1:=dt*y2;

dy2:=dt*ff(m,km,r,c,s,p,y1);

y1:=y1+dy1;

y2:=y2+dy2;

t:=t+dt;

end;

c:=c+0.2;

k:=k+1;

end;

//2

k:=0;

t:=30;y1:=vo;y2:=ff(m,km,r,c,s,p,vo);

while t<=44 do

begin

c:=0.1;

while c<=2 do

begin

form5.Chart2.SeriesList[k].Add(y2,floattostr(c));

dy1:=dt*y2;

dy2:=dt*ff(m,km,r,c,s,p,y1);

y1:=y1+dy1;

y2:=y2+dy2;

c:=c+0.1;

end;

t:=t+2;

k:=k+1;

end;

end;

//Кнопка ввода данных

procedure TForm4.N1Click(Sender: TObject);

//Кнопка построения графика

procedure TForm4.N2Click(Sender: TObject);

begin

//очистка старых графиков

for ind:=0 to 7 do

begin

form4.Chart1.SeriesList[ind].Clear;

form4.Chart2.SeriesList[ind].Clear;

end;

//Входные параметры

r:=strtofloat(edit1.Text);//радиус падующего тела

km:=strtofloat(edit2.text);//коэффициент вязкости

s:=strtofloat(edit3.Text);//площадь падующего тела

p:=strtofloat(edit4.Text);//плотность среды

ho:=strtofloat(edit5.Text);//начальная высота

vo:=strtofloat(edit6.Text);//начальная скорость

m:=strtofloat(edit7.Text);//масса падающего тела

cn:=0.1;ck:=1.5;//интервал лобового сопротивления среды

dc:=(1.5-0.1)/7;//шаг изменения лобового сопротивления

k:=0;//номер графика

// цикл по лобовому сопротивлению

while cn<=ck do

begin

v:=vo;//первая производная расстояния по времени- начальная скорость

a:=1; //вторая производная расстояния по времени-ускорение

ind:=1;//счетчик записи данных в массивы

t:=0;//счетчик подчета времени падения тела

h1:=ho;//высота с которой падает тело

h2:=0; // тело лежит на земле

//цикл времени падения тела

while h1>=h2 do

begin

//запись в массивы полученных данных

arrh[ind]:=h1;//высота

arrv[ind]:=v;//скорость

arrt[ind]:=t;//время

//вывод графиков

form4.Chart1.SeriesList[k].Add(arrh[ind],floattostr(arrt[ind]));

form4.Chart2.SeriesList[k].Add(arrv[ind],floattostr(arrt[ind]));

// решение дифференциального уравнения методом Эйлера

dv:=0.01*a;

da:=0.01*ff(m,km,r,cn,s,p,a);

v:=v+dv;

a:=a+da;

h1:=h1-arrv[ind]*arrt[ind];//изменение высоты за время t

t:=t+0.01;//изменение времени

ind:=ind+1;

end;

//переход к следующему коэффициенту лобового сопротивления

// и соответсвенно к следующему графику

k:=k+1;

cn:=cn+dc;

end;

end;

procedure TForm4.N4Click(Sender: TObject);

begin

form4.Hide;

fhello_html_m3952f3dc.pngorm5.hide;

form1.Show;

end;









hello_html_m1898225c.png
















hello_html_m112f8e22.png






















Пример2 . Исследовать падение шарика радиуса R с высоты h в среде MatLab

Задачи исследования:

1. Получить зависимости h=f(t), v=f(t) при начальных значениях r,v0,h0,k1,k2.

2. Получить семейство графиков h=f(t,r), v=f(t,r).

3. Получить зависимость h=f(r), v=f(r) при t - const

4. Исследовать в среде MatLab


hello_html_9c30296.gif





Исследование системы с использованием программы MatLab.

Графики зависимостей v(t), h(t) для r=5mm c

помощью решателей ode15,ode45


hello_html_62891098.gif



m-фа M-file f1.m:

function F1=pli(t,y)

mu=0.018; //динамическая вязкость воздуха

p1=800; //плотность шарика

p2=1.29; //плотность воздуха

c=0.4; //лобовое сопротивление шарика

g=9.8;

k=1000; //число разбиений

r=0.005; //радиус шарика

F1=[-y(2);9.8-9*mu/(r*r*p1*2)*y(2)-3*c*y(2)*y(2)*p2/(8*r*p1)];

m-файл для решения ДУ с помощью программы на языке программы MatLab M-file Sambo.m:

mu=0.018; //динамическая вязкость воздуха

p1=800; //плотность шарика

p2=1.29; //плотность воздуха

c=0.4; //коэффициент лобового сопротивления для шарика

g=9.8; //ускорение свободного падения

k=1000; //число разбиений

y(2)=0;y(1)=10;//начальные условия

x0=0;xk=5;//границы изменения времени

dx=(xk-x0)/k; //шаг интегрирования

r=0.005; //радиус шарика

x=x0;

subplot(2,1,1); //деление формы на 2 граф. окна

hold on; //обеспечивает продолжение вывода графиков в окно

for i=0:k-1

if y(1)<0 break,end;йл для решения ДУ с помощью решетелей

// вычисление коэффициентов для метода Рунге - Кутта

k11=-dx*y(2);k21=dx*(9.8-3*mu/(r*r*p1*2)-c*pi*y(2)*y(2)*p2/(8*r*p1));

k12=-dx*(y(2)+k21/2);

k22=dx*(9.8-3*mu/(r*r*p1*2)-c*pi*(y(2)+k21/2)*(y(2)+k21/2)*p2/(8*r*p1)); k13=-dx*(y(2)+k22/2); k23=dx*(9.8-3*mu/(r*r*p1*2)-c*pi*(y(2)+k22/2)*(y(2)+k22/2)*p2/(8*r*p1));

k14=-dx*(y(2)+k23);k24=dx*(9.8-3*mu/(r*r*p1*2)-c*pi*(y(2)+k23)

*(y(2)+k23)*p2/(8*r*p1)); //расчёт разности

d1=(k11+2*k12+2*k13+k14)/6; d2=(k21+2*k22+2*k23+k24)/6;

//вычисление новых значений у1 и у2

y(11)=y(1)+d1; y(21)=y(2)+d2;

plot([x x+dx],[y(1) y(11)],’r’);plot([x x+dx],[y(2) y(21)],’b’),

x=x+dx; //новое значение х

y(1)=y(11); y(2)=y(21);

end;

grid on; //добавляет сетку к текущему графику

title(“Kpuvue v(t) u h(t) npu r = 5mm”); //установка на

графике титульной надписи

legend(“vucota”,’ckopoct’); //добавление к текущему графику xlabel('t');ylabel('v,h');

plot([0 5],[0 0],'k');axis tight;hold off; y(2)=0;y(1)=10;x0=0;xk=10;dx=(xk-x0)/k;

r0=0.002;rk=0.005;dr=(rk-r0)/3;

n=2;x=x0;r=r0;

subplot(2,1,2);hold on;

while r<=rk

for i=0:k-1

if y(1)<0 break,end;

k11=-dx*y(2);

k21=dx*(9.8-3*mu/(r*r*p1*2)-c*pi*y(2)*y(2)*p2/(8*r*p1));

k12=-dx*(y(2)+k21/2);

k22=dx*(9.8-3*mu/(r*r*p1*2)-c*pi*(y(2)+k21/2)*(y(2)+k21/2)*p2/(8*r*p1)); k13=-dx*(y(2)+k22/2);

k23=dx*(9.8-3*mu/(r*r*p1*2)-c*pi*(y(2)+k22/2)*(y(2)+k22/2)*p2/(8*r*p1)); k14=-dx*(y(2)+k23);

k24=dx*(9.8-3*mu/(r*r*p1*2)-c*pi*(y(2)+k23)*(y(2)+k23)*p2/(8*r*p1));

d1=(k11+2*k12+2*k13+k14)/6; d2=(k21+2*k22+2*k23+k24)/6;

y(11)=y(1)+d1; y(21)=y(2)+d2;

plot([x x+dx],[y(1) y(11)],'r'); plot([x x+dx],[y(2) y(21)],'b');

x=x+dx; y(1)=y(11); y(2)=y(21); end;

x=x0; y(2)=0; y(1)=10;

r=r+dr; end;grid on;

title('Cemeuctvo v(t) u h(t): 2mm < r < 5mm');

legend('vucota','ckopoct');xlabel('t');ylabel('v,h');

plot([0 5],[0 0],'k');axis tight;hold off;

Графики зависимостей v(t),h(t) при r=5 mm и

семейство графиков для 2mm


hello_html_m2c84c8c2.gif



















hello_html_m7cf85562.gif

Simulink


hello_html_m5373e497.gif






Пример3. Исследовать движение исследовательского зонда вертикально вверх с летящего самолета

Условие

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

Входные параметры.

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


hello_html_6b5b4866.gif


Входные параметры модели:

¦ начальная высота тела;

¦ начальная скорость тела;

¦ величины, определяющие коэффициенты сопротивления среды

function FmO=fun(t,y)

tv=10;

a=100;

g=9,8

mn=10000;

mk=3000;

m=mn-a*t;

ifm>mk

else

m=mk end; mm=m; if t>=tv

p=0;

m=500; else

p=1000000;

m=mam; end;

FmO=[(p*cos(y(2))-cos(y(2))-g*sin(y(2)))/m;

((p*sin(y(2))+sin(y(2)))/m-

g*cos(y(2)))/y(l);-y(l)*cos(y(2))/10;y(l)*sin(y(2))/10];

Y0=[1000 pi/2 0 600];

[T,Y]=odel5s(@fun,[0 1000],Y0);

plot(Y(:53),Y(:,4),’g’);

hold on;

title(“zavisimost visoti ot vremeni”);

xlabel(fts’);

у1аЬеl(Ът’);

axis([0,5000,0,2000]);

hello_html_412de58b.gif








Пример. Исследование динамики объектов, брошенных под углом к горизонту.

Цель занятия:

Получить практические навыки исследования объектов, динамика которых описывается дифференциальными уравнениями 2-го порядка с использованием языков программирования :Delphi, VC++, VC# и информационных технологий MatLab.

Задачи занятия:

1. Разработка алгоритмов

2. Построение семейства кривых y(x), dy/dx(x) при параметрах a-const, a-var

3. Получение зависимостей y(a), dy/da при x-const

4. Анализ результатов исследований.

Задания для проведения лабораторной работы


Дhello_html_6b5b4866.gifинамика объекта описывается системой дифференциальных уравнений







mdvx/dt = -Fccosq; =-Fcvx/v0

mdvy/dt = -m*g - Fcsinq,=-m*g - Fcvy/v0,


где Fc =k1V+ k2*V2 - линейная и квадратичная составляющая сопротивления воздуха.

2. Динамика объекта может представляется и алгебраическими уравнениями вида

hello_html_200c990e.gif


hello_html_m2d8c2430.gif




var

Form1: TForm1;

xdp,ydp:integer;

x0,y0:integer;

alfa,veloc:real;

intravel:boolean=false;

// выводит результатскорости и направления

function rtostr(t:real):string;

var e:string;

begin

str(t:4:1,e);

result:=e;

end;

//выход

procedure TForm1.Button1Click(Sender: TObject);

begin

halt;

end;

// рисует стрелку , выводит скорость и направление

procedure drawmousec;

begin

with form1.paintbox1.Canvas do begin

pen.Mode:=pmXor;

moveto(x0,y0);

Lineto(xdp,ydp);

alfa:=pi/2-arctan((ydp-y0)/(xdp-x0));

Lineto(trunc(xdp+15*cos(-pi/2-alfa+pi/10)),trunc(ydp+15*sin(-pi/2-alfa+pi/10)));

moveto(xdp,ydp);

Lineto(trunc(xdp+15*cos(-pi/2-alfa-pi/10)),trunc(ydp+15*sin(-pi/2-alfa-pi/10)));

end;

alfa:=arctan((y0-ydp)/(xdp-x0));

form1.label1.Caption:='Направление: '+rtostr(alfa/pi*180)+' ; скорость: '+rtostr(veloc/50)+' ';

veloc:=sqrt(sqr(ydp/1-y0)+sqr(xdp/1-x0));

end;

// задает камень

procedure TForm1.PaintBox1Paint(Sender: TObject);

begin

if intravel then exit;

paintbox1.Canvas.pen.color:=clwhite;

paintbox1.Canvas.brush.color:=clwhite;

drawmousec;

end;

//координаты замка

procedure TForm1.FormCreate(Sender: TObject);

begin

xdp:=shape31.Left;

ydp:=shape31.Top;

x0:=shape25.left+shape25.width div 2;

y0:=shape25.top+shape25.height div 2;

end;

var dr:boolean=false;

//передвижение мышки вниз

procedure TForm1.PaintBox1MouseDown(Sender: TObject; Button: TMouseButton;

Shift: TShiftState; X, Y: Integer);

begin

if (x>x0)and(y

dr:=true;

drawmousec;

ydp:=y;xdp:=x;

drawmousec;

end;

end;

//передвижение мыши (стелки)

procedure TForm1.PaintBox1MouseMove(Sender: TObject; Shift: TShiftState; X,

Y: Integer);

begin

if (dr)and(x>x0)and(y

drawmousec;

ydp:=y;xdp:=x;

drawmousec;

end;

end;

// бросок

procedure TForm1.Button2Click(Sender: TObject);

const dt=0.01;

g=9.81;

var x,y,vx,vy,xo,yo:real;

bad:boolean;

r:integer;

begin

intravel:=true;

drawmousec;

vx:=veloc*cos(alfa);

vy:=veloc*sin(alfa);

x:=shape25.Left+shape25.height div 2;y:=shape25.Top+shape25.width div 2;r:=0;

repeat

xo:=x;yo:=y;

x:=x+vx*dt;y:=y-vy*dt;inc(r);

if r=10 then begin

sleep(1);

r:=0;

end;

vy:=vy-g*dt;

paintbox1.Canvas.pen.Mode:=pmCopy;

paintbox1.Canvas.pen.color:=clgreen;

paintbox1.Canvas.brush.color:=clgreen;

paintbox1.Canvas.Ellipse(trunc(xo-4),trunc(yo-4),trunc(xo+5),trunc(yo+5));

paintbox1.Canvas.pen.color:=clwhite;

paintbox1.Canvas.brush.color:=clwhite;

paintbox1.Canvas.Ellipse(trunc(x-3),trunc(y-3),trunc(x+4),trunc(y+4));

//trying to end the journey...

bad:=true;

if (y>y0/1)or(y<=0) then break;

if x>=shape15.Left then begin

if (x<=shape15.Left+shape15.width)and((y

(y>shape16.Top+shape16.height)) then break;

end;

if x>=shape22.Left then begin

if (y

(y>shape22.Top+shape22.height) then break else begin

bad:=false;

break;

end;

end;

until false;

intravel:=false;

if bad then begin

sleep(1000);

paintbox1.Repaint;

// drawmousec;

end else begin

shape32.Visible:=true;

label2.Visible:=true;

sleep(1000);

paintbox1.Repaint;

// drawmousec;

end;

end;

// изменение размеров замка

procedure TForm1.Shape32MouseDown(Sender: TObject; Button: TMouseButton;

Shift: TShiftState; X, Y: Integer);

begin

shape32.Visible:=false;

label2.Visible:=false;

end;

//вывод результатов

procedure TForm1.Label2MouseDown(Sender: TObject; Button: TMouseButton;

Shift: TShiftState; X, Y: Integer);

begin

shape32.Visible:=false;

label2.Visible:=false;

end;

// изменение координат замка

procedure TForm1.Splitter1CanResize(Sender: TObject; var NewSize: Integer;

var Accept: Boolean);

begin

accept:=false;

if newsize>=50 then begin

accept:=true;

end;

end;

// изменение координат замка

procedure TForm1.Splitter1Moved(Sender: TObject);

begin

halt;

end;

var s1d:boolean=false;

var s2d:boolean=false;

var s3d:boolean=false;

var s4d:boolean=false;

var s5d:boolean=false;

procedure TForm1.PaintBox2MouseDown(Sender: TObject; Button: TMouseButton;

Shift: TShiftState; X, Y: Integer);

begin

s1d:=true;

end;

// изменение координатов замка мышью

procedure TForm1.PaintBox2MouseMove(Sender: TObject; Shift: TShiftState; X,

Y: Integer);

begin

if (s1d)and(shape15.Left+x>=50)and(shape15.Left+shape15.Width+x+50<=shape21.left) then begin

paintbox2.;

paintbox3.;

paintbox4.;

shape15.;

shape16.;

shape17.;

shape18.;

shape19.;

shape20.;

shape23.;

shape23.Width:=form1.width;

end;

end;

//передвижение мыши для изменения координат замка вверх

procedure TForm1.PaintBox2MouseUp(Sender: TObject; Button: TMouseButton;

Shift: TShiftState; X, Y: Integer);

begin

s1d:=false;

end;

//передвижение мыши для изменения координат замка вниз

procedure TForm1.PaintBox3MouseDown(Sender: TObject; Button: TMouseButton;

Shift: TShiftState; X, Y: Integer);

begin

s2d:=true;

end;

// изменение координатов замка мыщью

procedure TForm1.PaintBox3MouseMove(Sender: TObject; Shift: TShiftState; X,

Y: Integer);

begin

if (s2d)and(shape15.Width+x>=10)and(shape15.Left+shape15.Width+x+50<=shape21.left) then begin

shape15.Width:=shape15.Width+x;

shape16.Width:=shape16.Width+x;

paintbox3.;

shape17.width:=shape15.width+20;

shape17.;

paintbox4.width:=shape15.width+20;

paintbox4.;

shape18.width:=shape17.width div 5;

shape19.width:=shape17.width div 5;

shape18.;

shape19.;

shape20.;

shape20.width:=shape17.left+shape17.width-shape20.left;

end;

end;

//передвижение мыши для изменения координат замка вверх

procedure TForm1.PaintBox3MouseUp(Sender: TObject; Button: TMouseButton;

Shift: TShiftState; X, Y: Integer);

begin

s2d:=false;

end;

//передвижение мыши для изменения координат замка вниз

procedure TForm1.PaintBox4MouseDown(Sender: TObject; Button: TMouseButton;

Shift: TShiftState; X, Y: Integer);

begin

s3d:=true;

end;

//передвижение мыши для изменения координат замка вверх

procedure TForm1.PaintBox4MouseUp(Sender: TObject; Button: TMouseButton;

Shift: TShiftState; X, Y: Integer);

begin

s3d:=false;

end;

//передвижение мыши для изменения координат

procedure TForm1.PaintBox4MouseMove(Sender: TObject; Shift: TShiftState; X,

Y: Integer);

begin

if (s3d)and(shape17.top+y>=100)and(shape16.Top+shape16.Height+y+20<=shape23.Top) then begin

shape17.;

shape18.;

shape19.;

shape16.;

shape20.;

shape15.;

paintbox4.;

paintbox3.;

paintbox2.;

shape15.height:=form1.height-shape15.top;

paintbox2.height:=form1.height-paintbox2.top;

paintbox3.height:=form1.height-paintbox3.top;

end;

end;

//передвижение мыши для изменения координат замка вниз

procedure TForm1.PaintBox5MouseDown(Sender: TObject; Button: TMouseButton;

Shift: TShiftState; X, Y: Integer);

begin

s4d:=true;

end;

//передвижение мыши для изменения координат замка вверх

procedure TForm1.PaintBox5MouseUp(Sender: TObject; Button: TMouseButton;

Shift: TShiftState; X, Y: Integer);

begin

s4d:=false;

end;

//передвижение мыши для изменения координат замка

procedure TForm1.PaintBox5MouseMove(Sender: TObject; Shift: TShiftState; X,

Y: Integer);

begin

if (s4d)and(shape22.top+y-10>=shape21.top)and(shape22.Top+shape22.Height+y+10<=shape23.Top) then begin

shape22.;

paintbox5.;

paintbox6.;

end;

end;

//передвижение мыши для изменения координат замка вниз

procedure TForm1.PaintBox6MouseDown(Sender: TObject; Button: TMouseButton;

Shift: TShiftState; X, Y: Integer);

begin

s5d:=true;

end;

//передвижение мыши для изменения координат замка вверх

procedure TForm1.PaintBox6MouseUp(Sender: TObject; Button: TMouseButton;

Shift: TShiftState; X, Y: Integer);

begin

s5d:=false;

end;

procedure TForm1.PaintBox6MouseMove(Sender: TObject; Shift: TShiftState; X,

Y: Integer);

begin

if (s5d)and(shape22.height+y>=25)and(shape22.Top+shape22.Height+y+10<=shape23.Top) then begin

shape22.height:=shape22.height+y;

paintbox6.;

end;

end;


procedure TForm1.PaintBox1MouseUp(Sender: TObject; Button: TMouseButton;

Shift: TShiftState; X, Y: Integer);

begin

dr:=false;

end;


end.



hello_html_2fffce59.png


















hello_html_11217e68.png



















Пример модели движения небесных тел

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

hello_html_m538e7452.gif



В этой формуле G=6,67*10-11 m3/(kг*с2) -гравитационная постоянная. В данной системе координат начало координат расположено на одном теле M. Уравнения, описывающие движение тела m в указанной системе координат, имеет вид

hello_html_m657fe4ef.gif



В проекциях на оси координат уравнение можно привести с следующей системе из 4-х дифференциальных уравнений



hello_html_m742e4bcd.gif


















Пhello_html_m4f243884.gifример 2. Исследовать выполнение второго закона Кеплера при движении различных планет.



hello_html_50da9bbb.gif


















unit Unit2;

var

Form2: TForm2;

angle,j,k,kx,ky,ra,rb,e,x,y,va,vb,vc,a,t,t2,speed:real;

xc,yc,w,h,xf,n,cc,step:integer; flag:boolean;

implementation

uses Unit1, Unit3;

{$R *.dfm}

procedure TForm2.FormClose(Sender: TObject; var Action: TCloseAction);

begin

form1.close;

end;

procedure TForm2.BitBtn1Click(Sender: TObject);


var rrr:trect; aa,bb,cc:integer; s:string;

begin

// дано

val(form1.edit1.text,ra,cc);

if cc<>0 then ra:=1; val(form1.edit2.text,e,cc);

if cc<>0 then e:=0.0167; val(form1.edit3.text,vc,cc);

if cc<>0 then vc:=0.9856;

t:=360/vc; t2:=t/2; flag:=true;

// вычисляем

va:=vc*sqrt((1+e)/(1-e)); vb:=vc*sqrt((1-e)/(1+e));

a:=(va-vb)/t2; val(edit1.Text,step,cc);

if cc<>0 then step:=1; str(va:3:11,s);

label1.caption:=s; str(vb:3:11,s);

label2.caption:=s; str(vc:3:11,s);

label3.caption:=s; rb:=ra*sqrt(1-e*e);

// графика

w:=image1.Width; h:=image1.Height;

xc:=w div 2; yc:=h div 2;

kx:=ra/(xc-50); ky:=rb/(yc-50);

if kx>ky then k:=(xc-20)/ra else k:=(yc-20)/rb;

timer1.Enabled:=true;

aa:=round(ra*k); bb:=round(rb*k);

image1.Canvas.Brush.Color:=clnavy;

rrr:=rect(0,0,w,h); image1.Canvas.FillRect(rrr);

image1.Canvas.Brush.Color:=clwhite;

rrr:=rect(xc-aa,yc-bb,xc+aa,yc+bb);

image1.Canvas.Ellipse(rrr);

cc:=3; image1.Canvas.Brush.Color:=clnavy;

rrr:=rect(xc-aa+cc,yc-bb+cc,xc+aa-cc,yc+bb-cc);

image1.Canvas.Ellipse(rrr); xf:=round(ra*e*k);

cc:=30; image1.Canvas.Brush.Color:=clyellow;

rrr:=rect(xc+xf+cc,yc+cc,xc+xf-cc,yc-cc);

image1.Canvas.Ellipse(rrr);

cc:=6; image1.Canvas.Brush.Color:=clred;

rrr:=rect(xc-xf+cc,yc+cc,xc-xf-cc,yc-cc);

image1.Canvas.Ellipse(rrr);

rrr:=rect(xc+xf+cc,yc+cc,xc+xf-cc,yc-cc);

image1.Canvas.Ellipse(rrr);

angle:=0; n:=0; speed:=va; j:=pi/180;

end;

procedure TForm2.Timer1Timer(Sender: TObject);

var rrr:trect;s:string;i:integer;

begin

image1.Canvas.MoveTo(xc+xf,yc);

image1.Canvas.pen.Color:=clyellow;

ifangle<360 then image1.Canvas.

LineTo (xc+round(ra*cos(angle*j)*k),

(c-round(rb*sin(angle*j)*k)); for i:=1 to step do

begin

inc(n);

if angle<180 then speed:=speed-a else speed:=speed+a;

angle:=angle+speed; end;

str(speed:3:11,s); label4.Caption:=s;

if angle>=(360-a) then begin timer1.Enabled:=false; flag:=false; end;

cc:=6; image1.Canvas.Brush.Color:=clred;

rrr:=rect(xc-xf+cc,yc+cc,xc-xf-cc,yc-cc);

image1.Canvas.Ellipse(rrr);

rrr:=rect(xc+xf+cc,yc+cc,xc+xf-cc,yc-cc);

image1.Canvas.Ellipse(rrr);end;

procedure TForm2.BitBtn2Click(Sender: TObject);

begin

if flag then timer1.enabled:=not timer1.enabled;

end;

procedure TForm2.BitBtn3Click(Sender: TObject);

begin

form2.Hide; form1.Show;end;

procedure TForm2.BitBtn4Click(Sender: TObject);

begin

form2.Hide; form3.Show;

end;end.


hello_html_m8375deb.gif




2.4 Лабораторные работы по разработке имитационных моделей


Пример. Разработка информационной модели студента ( учащихся)

Цель занятия: Получить практические навыки в разработке информационных систем в среду БД Access , Excel

Задачи занятия:

2. Создать графическую модель студента с отображением диаграммы оценок по темам конкретной дисциплины

Информационная модель студента

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

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

-стабильность обучения по отдельным предметам, темам и дисциплинам;

- прогнозирование результатов обучения в предстоящем семестре;

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

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

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


Пример имитационной модели знаний студентов в среде Excel









hello_html_m61a6a2c6.png











hello_html_673cb3fa.png











hello_html_44c41c09.png

hello_html_m3d678ce0.pnghello_html_37771e0f.png


















hello_html_4bf03a05.png


hello_html_m22869212.png


hello_html_5aed437e.png











hello_html_m39abbf24.png






hello_html_m20c94188.png












hello_html_m49f3c373.png















Пример имитационной модели учащихся школы в среде Access

hello_html_m38fb1c54.gif








hello_html_m4daf4484.gif



hello_html_2b316123.gif



hello_html_4188bf1c.gif












hello_html_1a032e94.gif













hello_html_5d061bbe.gif













hello_html_31833f86.gif












hello_html_62a648b3.gif














hello_html_mb387cf6.gif












hello_html_7374c5b9.gif













2.5 Разработка моделей транспортных задач

Пример. «Размещения предприятий»


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

1.(перевозится неотрицательное количество продукта);

2. (выпускаемое количество продукта не

больше возможного обьема производства и равно вывозимому количеству продукта);

3. , j=1,...,n (потребности потребителей),

(каждый пункт потребления получает столько, сколько ему требуется).

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

hello_html_18870d0f.gifhello_html_m2df761c1.gifhello_html_m2df761c1.gifhello_html_mfd55774.gif



hello_html_5cafc885.gif

,

где - затраты на производство,

- затраты на транспортировку.

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

,

где - затраты на производство,

- затраты на транспортировку.

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








Пример Моделирование системы планирования

на основе метода сетевого графа

Задачи исследования

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

Исходные данные



Задание: Дана схема, нужно найти критический путь



hello_html_7f8d5090.png






hello_html_m131adbf8.png




hello_html_m38651be6.png

















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

Условие. Пусть предприятие производит столы и стулья. Расход ресурсов на их производство и прибыль от их реализации представлены в таблице. столы стулья объем ресурсов

Расход древесины

на изделие, м3 0,5 0,04 200

Расход труда(чел-час) 12 0,6 1800


Прибыль от реализации

единицы изделия , руб 180 20

Кроме того на производство 80 столов имеется заказ с министерством

Уравнения1-количество столов х2 - количество стульев

0,5 х1 + 0,04х2<=200

12x1 + 0,6x2 <=18000

x1+20x2 =max ( целевая функция)


hello_html_m775891b7.png




















hello_html_63312721.png

















hello_html_1c242426.png









hello_html_59ef985b.gif

















hello_html_4bc1a142.gif












hello_html_29ac9ad5.png




unit











hello_html_m4e3a6278.png























Unit4;

procedure TForm4.Button3Click(Sender: TObject);

begin

close();

end;

procedure TForm4.Button1Click(Sender: TObject);

begin

form1.visible:=true;

form4.visible:=false;

end;

function Fy(x:real): real;

begin

Fy := (60-6*x)/4;

end;

function Fx(x:real): real;

begin

Fx :=(100+2*x)/10;

end;

procedure TForm4.Button2Click(Sender: TObject);

var a, b ,h, x1, y, x2: Double;

begin

Chart1.Series[0].Clear;

Chart1.Series[1].Clear;

Chart1.Series[2].Clear;

a := StrToFloat(LabeledEdit1.Text); //левая граница

b := StrToFloat(LabeledEdit2.Text); //правая граница

h := StrToFloat(LabeledEdit3.Text); //шаг

x1 := a;

repeat

y := Fy(x1);

x2 := Fx(y);

hart1.Series[0].AddXY(x1,y,FloatToStrF(x1,ffGeneral,4,3),clRed); //график F(x)

hart1.Series[1].AddXY(x2,y,FloatToStrF(y,ffGeneral,4,3),clGreen); //график F(y)

x1 := x1 + h; //приращение координат

until x1 > b;

x1 := a;

repeat

y := Fy(x1);

x2 := Fx(y);

x1 := x1 + h/10000;

if ( x1 < x2 + h/10000 ) and ( x1 > x2 - h/10000 ) then

begin

Label5.Caption := 'Решение: ' +#13+ 'x = '+FloatToStrF(x1,ffGeneral,7,7)+#13;

Label5.Caption :=Label5.Caption + 'y = '+FloatToStrF(y,ffGeneral,7,7);

Chart1.Series[2].AddXY(x1,y,FloatToStrF(x1,ffGeneral,4,3),clBlue);

end;

until x1 > b;

end;end.





2.9 Лабораторная работа

Тема: Исследование вероятностных систем ( систем массового обслуживания) с использованием различных информационных технологий.

Цель занятия: Получить практические навыки исследования систем ( объектов) с вычислением параметров статистических моделей и определением эффективности исследуемых систем с использованием языков программирования : Delphi, VC++, VC# и информационных технологий MatLab, Scilab..

Задачи занятия:

1. Разработка алгоритмов

2. Разработка программы исследования различных систем

3. Анализ результатов исследований с вычислением параметров и эффективности систем.


Задания для исследований в разделе Задания

Пhello_html_m2840eacf.pngример 1. Исследовать параметры обслуживания покупателей в магазине с одним продавцом





hello_html_m58344ca6.gif





hello_html_m418a363d.png


Программа на Делфи

unit Unit2;

function f(x:real;g:real;h:real):real;

begin

f:=(1/(g*sqrt(3.14*2)))*(power(exp(1),(((x-h)*(x-h))/(-2*g*g)))); // нормальный закон распределения

end;

procedure TForm2.Button2Click(Sender: TObject);

var n,i,a,b,rk:integer;

g,h,f1,fn,h1,x1,sr,t,y,kol,vrop:real;

begin

randomize;

g:=strtofloat(edit1.text); // среднеквадратич отклонение X

h:=strtofloat(edit2.text);//мат ожидание

Chart1.Series[0].Clear; //очистка графика

a:=StrToint(Edit4.Text); //левая граница

b:=StrToint(Edit5.Text); //правая граница

h1:=0.1;

x1:=a;

n:=1;

fn:=0;

t:=h1;

stringgrid1.Cells[0,0]:='x';

stringgrid1.Cells[1,0]:='f(x)';

repeat

f1:=f(x1,g,h);

Chart1.Serieslist[0].AddXY(x1,f1,'',clRed); //график F(x)

stringgrid1.Cells[0,n]:=floattostr(x1);

stringgrid1.Cells[1,n]:=floattostr(f1);

x1:=x1+h1;

n:=n+1;

fn:=f1+fn;

stringgrid1.RowCount:=n+1;

until x1>b;

kol:=(b-a)*60*g;

edit6.text:=floattostr(kol);

rk:=strtoint(edit6.text);

i:=random(rk);

edit7.text:=inttostr(random(5));

sr:=strtofloat(edit7.text);

vrop:=((b-a)*60-sr*kol)/kol;

if vrop>=0 then

edit8.Text:=floattostr(((b-a)*60-sr*kol)/kol)

else edit8.Text:='0';

if vrop<=0 then

edit3.Text:=floattostr(-vrop) else

edit3.Text:='0';

edit6.text:=floattostr(i);

end;

end.

















hello_html_m3cbf2220.gif




























hello_html_40780ef1.png






2.9











hello_html_4f90fa18.png


















hello_html_m7767dd41.png

















unit Unit3;

interface

uses

Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms, Dialogs, Menus, ExtCtrls, StdCtrls, jpeg;

type

TForm3 = class(TForm)

MainMenu1: TMainMenu;

N1: TMenuItem;

N2: TMenuItem;

Timer1: TTimer;

GroupBox1: TGroupBox;

Edit1: TEdit;

Edit2: TEdit;

Edit3: TEdit;

Edit4: TEdit;

Edit5: TEdit;

Label1: TLabel;

Label2: TLabel;

Label3: TLabel;

Label4: TLabel;

Label5: TLabel;

Timer2: TTimer;

Image1: TImage;

N3: TMenuItem;

Image2: TImage;

Label6: TLabel;

procedure N2Click(Sender: TObject);

procedure Timer1Timer(Sender: TObject);

procedure N1Click(Sender: TObject);

procedure Timer2Timer(Sender: TObject);

procedure N3Click(Sender: TObject);

var

Form3: TForm3;

a,b,t,kol,za,k,i,l,j,ll,gdet,zanvr:Integer;

nach,vrem,kon,pokgdet,pok,x1,y1:Integer;

Buyx,Buyy:Integer;

obsl,prish:Integer;

tekvr,chas,minut:String;

pokin:array[1..1000] of real;

pokvr1:array[1..1000] of real;

pokvr:array[1..1000] of real;

prod:array[1..1000] of real;

ogid:array[1..1000] of real;

ob:array[1..1000] of real;

p,s,inte,buy,x,y,maxt,tt,vr:real;

implementation

uses Unit2, Unit1;

procedure TForm3.N2Click(Sender: TObject);

begin

form3.Hide;

form1.show;

end;

procedure TForm3.Timer1Timer(Sender: TObject);

begin

If zanvr=prod[k] Then {если пок-ля обслужили}

begin

timer2.Enabled:=true;

image1.Visible:=true;

timer2.Interval:=100;

zanvr:=0;

k:=k+1;

If ogid[k]=0 Then {если пок-ль не в очереди}

begin

za:=za+1; {увеличиваем кол-во пришедших}

edit4.Text:=inttostr(za);

end;

If kol>0 Then begin

kol:=kol-1; {уменьшаем кол-во ждущих в очереди}

edit5.Text:=inttostr(kol);

form3.Canvas.Pen.Color:=clskyblue;

form3.Canvas.Brush.Color:=clskyblue;

form3.Canvas.Ellipse(x1,y1,x1+30,y1+30);

x1:=x1-32;

end;

obsl:=obsl+1; {увел-м кол-во обслуженных}

edit2.Text:=inttostr(obsl);

end;

buy:=pokvr[k];

If buy=i Then begin {если покупатель пришел}

zanvr:=zanvr+1; {начинаем обслуживать}

pokvr[k]:=i+1;

end

Else gdet:=gdet+1; {иначе увел-ся простой продавца}

For j:=k+1 To l-1 do

begin

x:=pokvr[j];

y:=pokvr[k];

If x

pokvr[j]:=pokvr[k];

If ogid[j]=0 Then begin

kol:=kol+1; {увел-м кол-во ждущих}

edit5.Text:=inttostr(kol);

x1:=x1+32;

form3.Canvas.Pen.Color:=clmaroon;

form3.Canvas.Brush.Color:=clmaroon;

form3.Canvas.Ellipse(x1,y1,x1+30,y1+30);

{увел-м кол-во пришедших}

za:=za+1;

edit4.Text:=inttostr(za);

ogid[j]:=1;

end;

end;

end;

i:=i+1; {увеличиваем время}

{преобразуем время в часы и минуты}

str((i div 60),chas);

str((i mod 60),minut);

tekvr:=chas+':'+minut;

edit1.Text:=tekvr;

edit3.Text:=inttostr(gdet);

if i>=kon then timer1.Enabled:=false;

end;


procedure TForm3.N1Click(Sender: TObject);

var stroka,stroka1:string;

begin

For i:=1 To obsl do

begin

ob[i]:=pokvr[i]-prod[i]-pokvr1[i];

s:=s+ob[i];

end;

For i:=obsl+1 To l-1 do {суммируем время ожидания в оч-ди}

begin

ob[i]:=pokvr[i]-pokvr1[i];

s:=s+ob[i];

end;

s:=s/(l-1); {находим среднее арифметич-ое}

stroka1:=inttostr(gdet);

stroka:=inttostr(round(s));

showmessage('Среднее время ожидания в очереди покупателем равно'+stroka+'минут; простой продавца в ожидании прихода покупателей'+stroka1+'минут');

end;


procedure TForm3.Timer2Timer(Sender: TObject);

begin

Image1.;

If Image1.Left<112 Then begin

Timer2.Enabled:=False;

Image1.Visible:=False;

Image1.;

End;

end;


procedure TForm3.N3Click(Sender: TObject);

begin

nach:=strtoint(form2.edit1.Text)*60;

inte:=strtofloat(form2.Edit3.Text);

kon:=strtoint(form2.edit2.Text)*60;

pok:=0;

gdet:=1;

k:=0;

zanvr:=0;

vrem:=0;

obsl:=0;

i:=1;

vr:=nach;

gdet:=0;

k:=1;

x1:=290;

y1:=140;

maxt:=random(5); {максимальное время}

maxt:=(maxt+1)/100;

While vr

begin {интервал прихода покупателей}

pokin[i]:=random(7); {в соответствии с функцией Пуассона}

pokin[i]:=Int((1/inte)*((pokin[i]+1)/10))+1;

tt:=(inte*pokin[i])*(exp((-inte*pokin[i])*ln(2.71)));

If tt>maxt Then

begin

vr:=vr+pokin[i];

pokvr[i]:=vr; {время прихода пок-лей}

pokvr1[i]:=vr;

i:=i+1;

end;

end;

For l:=1 To i do

begin {время, кот. тратит продавец}

prod[l]:=random(8);

prod[l]:=Int((1/inte)*((prod[l]+1)/10))+1;

ob[l]:=0;

ogid[l]:=0;

end;

i:=nach;

timer1.enabled:=true;

timer1.interval:=500;

end;

end.
























Тема: Исследование оптимальных ( имитационных) систем Цель занятия: Получить практические навыки исследования оптимальных ( имитационных) систем

Задачи занятия:

1. Разработка алгоритмов

2. Разработка программы исследования различных систем

3. Анализ результатов исследований с вычислением параметров


hello_html_213a60a3.png





hello_html_m14cd8377.png







hello_html_8cf1de9.png













unit Unit1;

interface

usesWindows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms, Dialogs, StdCtrls;

var

Form1: TForm1;

dx,dy:real; //масштаб реки K_Der,D_Reki,N_R,D_Reki1,N_R1,R,i,num,mmx,rast,j,mmn,x,y,xv,yv,xn,yn:integer;

St:array[1..800,1..2] of integer; //длина реки

Pl:array[1..100,1..2] of integer; //количество деревень

Rw:array[1..2,1..2] of integer; //река

implementation

uses Unit2;

procedure TForm1.Button1Click(Sender: TObject);

begin

form1.Hide;

form2.show;

end;

procedure TForm1.Button2Click(Sender: TObject);

begin

xv:=strtoint(inputbox('vvod ','1-ая координата левого верхнего угла поля(0-10)',''));

yv:=strtoint(inputbox('vvod ','2-ая координата левого верхнего угла поля(0-10)',''));

xn:=strtoint(inputbox('vvod ','1-ая координата правого нижнего угла поля(580-650)',''));

yn:=strtoint(inputbox('vvod ','2-ая координата правого нижнего угла поля(500-580)',''));

K_Der:=strtoint(inputbox('vvod ','количество деревень(1-70)',''));

D_Reki:=strtoint(inputbox('vvod ','длина реки(300-550)',''));

N_R:=strtoint(inputbox('vvod ','начало реки(10-200)',''));

{D_Reki1:=strtoint(inputbox('vvod ','длина реки(300-550)',''));

N_R1:=strtoint(inputbox('vvod ','начало реки(10-200)',''));}

R:=strtoint(inputbox('vvod ','длина разбиения(10-20)',''));

randomize;

for i:=1 to K_Der do begin //определение координат деревень

Pl[i,1]:=random(570)+10;

Pl[i,2]:=random(480)+10;

end;

Rw[1,1]:=N_R; //координаты реки

Rw[1,2]:=D_Reki;

Rw[2,1]:=random(480)+10;

Rw[2,2]:=random(480)+10;

end;

procedure TForm1.Button3Click(Sender: TObject);

begin

form1.Canvas.Brush.Color:=clolive; //вывод поля

form1.Canvas.Rectangle(xv,yv,xn,yn);

dx:=D_Reki/(D_Reki-1); //определение масштаба по реке

dy:=(Rw[2,2]-Rw[2,1])/(D_Reki-1);

form1.Canvas.Pen.Color:=clblue;

form1.Canvas.MoveTo(Rw[1,1],Rw[2,1]-3); //вывод реки на экран

form1.Canvas.LineTo(Rw[1,2],Rw[2,2]-3);

form1.Canvas.MoveTo(Rw[1,1],Rw[2,1]+3);

form1.Canvas.LineTo(Rw[1,2],Rw[2,2]+3);

for j:=1 to D_Reki do begin //определение расстояния до деревень

mmx:=0;

for i:=1 to K_Der do begin

rast:=round(sqrt(sqr(10+ dx*(j-1)-Pl[i,1])+sqr(Rw[2,1]+dy*(j-1)-Pl[i,2])));

if mmx

mmx:=rast; //максимальное из расстояний до деревень

St[j,1]:=rast;

end;end;

end;

for i:=1 to K_Der do begin

form1.Canvas.Brush.Color:=clyellow; //вывод деревень с их нумерацией

form1.Canvas.Pen.Color:=clred;

form1.Canvas.TextOut(Pl[i,1]-14,Pl[i,2]-14,inttostr(i));

end;

end;

procedure TForm1.Button4Click(Sender: TObject);

begin

{mmn:=1000;

for i:=1 to (D_Reki-N_R) do begin}

while i<=(D_Reki-N_R) do begin

mmn:=1000;i:=i+R; //определение искомого места на реке

if St[i,1]

mmn:=St[i,1];

num:=i+R;

end;end;

for j:=1 to K_Der do begin //вывод расстояний от искомого места до деревень

rast:=round(sqrt(sqr(10+ dx*(num-1)-Pl[j,1])+sqr(Rw[2,1]+dy*(num-1)-Pl[j,2])));

form1.Canvas.Brush.Color:=clsilver;

form1.Canvas.Pen.Color:=clred;

form1.Canvas.TextOut(Pl[j,1]-3,Pl[j,2]-3,inttostr(rast));

end;

form1.Canvas.Pen.Color:=clfuchsia; //вывод искомого места

form1.Canvas.brush.Color:=clfuchsia;

form1.Canvas.Ellipse(round(dx*(num-R)+2),round(Rw[2,1]+dy*(num-R)-8),round(18+dx*(num-R)),round(Rw[2,1]+dy*(num-R)+8));

for i:=1 to K_Der do begin //линии от искомого места до деревень

form1.Canvas.Pen.Color:=clnavy;

form1.Canvas.brush.Color:=clnavy;

form1.Canvas.MoveTo(round(dx*(num-1)+7),round(Rw[2,1]+dy*(num-1)-5));

form1.Canvas.LineTo(Pl[i,1]-5,Pl[i,2]-5);end;

end;end.











hello_html_m32baf86d.png


hello_html_58dfac1d.png



















hello_html_m1e6242ca.png


hello_html_7dbe1014.png


unit Unit4;


interface


uses

Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,

Dialogs, StdCtrls, ExtCtrls;


type

TForm4 = class(TForm)

Button1: TButton;

GroupBox1: TGroupBox;

PaintBox1: TPaintBox;

LabeledEdit1: TLabeledEdit;

ListBox1: TListBox;

ListBox2: TListBox;

Button2: TButton;

Label1: TLabel;

Label2: TLabel;

GroupBox2: TGroupBox;

RadioButton1: TRadioButton;

RadioButton2: TRadioButton;

LabelMouse: TLabel;

procedure Button1Click(Sender: TObject);

procedure Button2Click(Sender: TObject);

procedure RadioButton2Click(Sender: TObject);

procedure RadioButton1Click(Sender: TObject);

procedure LabeledEdit1Change(Sender: TObject);

procedure PaintBox1MouseMove(Sender: TObject; Shift: TShiftState; X,

Y: Integer);

procedure PaintBox1MouseUp(Sender: TObject; Button: TMouseButton;

Shift: TShiftState; X, Y: Integer);

procedure FormMouseMove(Sender: TObject; Shift: TShiftState; X,

Y: Integer);

procedure FormShow(Sender: TObject);

private

{ Private declarations }

public

{ Public declarations }

end;


Type

PointType = record

x: Integer;

y: Integer;

end;


var

Form4: TForm4;


ReadMouse: Boolean;

MouseClickCount: Byte;

n: integer; { количество деревень }

i, j: integer; { для цикла }

MinOfAverages: integer; { минимальный средний путь до всех деревень }

iMinOfAverages: integer;

Average, AverVal: integer;

point: array[1..50] of PointType; { для записывания координат деревень }

z: array[1..50,1..50] of integer; { для записывания расстояний между деревнями }

Mx, My: integer;


implementation


uses Unit1;


{$R *.dfm}


procedure Vill( vX, vY, vI: Integer; clColor: TColor);

begin

Form4.PaintBox1.Canvas.Pen.Color := clColor;

Form4.PaintBox1.Canvas.Brush.Color := clColor;

Form4.PaintBox1.Canvas.Ellipse( vX - 10, vY - 10,

vX + 10, vY + 10 );

Form4.PaintBox1.Canvas.TextOut( vX-5, vY-6, IntToStr(vI) );

Form4.PaintBox1.Canvas.Brush.Color := clBtnFace;

Form4.PaintBox1.Canvas.TextOut( vX - 15, vY + 11,

IntToStr(vX) + ',' + IntToStr(vY) );

end;


//кнопка "Выполнить"

procedure TForm4.Button1Click(Sender: TObject);

begin

Form4.PaintBox1.Repaint;

if length( LabeledEdit1.Text ) > 0 then

n := StrToInt( LabeledEdit1.Text );

ListBox2.Clear;

{ цикл вычисления растояний между деревнями }

for i := 1 to n - 1 do

begin

for j := i + 1 to n do

begin

z[i,j] := trunc( sqrt( sqr( point[i].x - point[j].x ) +

sqr( point[i].y - point[j].y ) ) );

z[j,i] := z[i,j];

ListBox2.Items.Add( 'Путь между деревнями ' + IntToStr(i) + ' и ' + IntToStr(j) + ' : ' + IntToStr(z[i,j]) + ' км.' );

end;

z[i,i] := 0;

end;


MinOfAverages := 10000;

ListBox1.Clear;

{ вычисление пункта,из которого средний путь до всех деревень минимален }

for i := 1 to n do

begin

PaintBox1.Canvas.Pen.Color := clBlue;

AverVal := 0;

for j := 1 to n do

begin

AverVal := AverVal + z[i,j];

PaintBox1.Canvas.MoveTo( point[i].x, point[i].y );

PaintBox1.Canvas.LineTo( point[j].x, point[j].y );

end;

Average := trunc( AverVal / ( n - 1 ) );

PaintBox1.Canvas.Pen.Color := clYellow;

ListBox1.Items.Add( 'Среднй путь от деревни № ' + IntToStr(i) + ' : ' + IntToStr(Average) + ' км.' );

if Average < MinOfAverages then

begin

MinOfAverages := Average;

iMinOfAverages := i;

end;

end;

Label1.Caption := 'Минимальный средний путь: ' + IntToStr(MinOfAverages);

Label2.Caption := 'Номер деревни,где будет располагаться больница: ' + IntToStr(iMinOfAverages);

PaintBox1.Canvas.Pen.Color := clGreen;

PaintBox1.Canvas.Brush.Color := clGreen;


for i := 1 to n do

Vill( point[i].x, point[i].y, i, clGreen);


{ Найденный пункт выделяем другим цветом }

Vill( point[iMinOfAverages].x, point[iMinOfAverages].y, iMinOfAverages, clRed );

{ PaintBox1.Canvas.Pen.Color := clRed;

PaintBox1.Canvas.Brush.Color := clRed;

PaintBox1.Canvas.Ellipse( point[iMinOfAverages].x-10,

point[iMinOfAverages].y-10,

point[iMinOfAverages].x+10,

point[iMinOfAverages].y+10 );

PaintBox1.Canvas.TextOut( point[iMinOfAverages].x - 5,

point[iMinOfAverages].y - 6,

IntToStr( iMinOfAverages ) ); }

end;

//кнопка "Закрыть"

procedure TForm4.Button2Click(Sender: TObject);

begin

form4.Hide;

form1.show;

end;


procedure TForm4.RadioButton2Click(Sender: TObject);

begin

PaintBox1.Repaint;

MouseClickCount := 0;

ReadMouse := True;

Button1.Enabled := False;

end;


procedure TForm4.RadioButton1Click(Sender: TObject);

begin

PaintBox1.Repaint;

ReadMouse := False;

Button1.Enabled := True;


Randomize;

{ цикл генерирования и вывода координат деревень }

for i := 1 to n do

begin

point[i].x := trunc( random( Mx ) );

point[i].y := trunc( random( My ) );

Vill( point[i].x, point[i].y, i, clGreen );

end;


end;


procedure TForm4.LabeledEdit1Change(Sender: TObject);

begin

if length( LabeledEdit1.Text ) > 0 then

n := StrToInt( LabeledEdit1.Text );

PaintBox1.Repaint;

if RadioButton1.Checked then RadioButton1Click(Sender)

else RadioButton2Click(Sender);


end;


procedure TForm4.PaintBox1MouseMove(Sender: TObject; Shift: TShiftState; X,

Y: Integer);

begin

if ReadMouse then LabelMouse.Visible := True;

LabelMouse.Caption := IntToStr(X) + ',' + IntToStr(Y);

end;


procedure TForm4.PaintBox1MouseUp(Sender: TObject; Button: TMouseButton;

Shift: TShiftState; X, Y: Integer);

begin

if ReadMouse then

begin

MouseClickCount := MouseClickCount + 1;

point[MouseClickCount].x := X;

point[MouseClickCount].y := Y;

Vill( X, Y, MouseClickCount, clGreen );

if MouseClickCount >= n then

begin

ReadMouse := False;

Button1.Enabled := True;

end;

end;

end;


procedure TForm4.FormMouseMove(Sender: TObject; Shift: TShiftState; X,

Y: Integer);

begin

LabelMouse.Visible := False;

end;


procedure TForm4.FormShow(Sender: TObject);

begin

Mx := PaintBox1.Width;

My := PaintBox1.Height;

Label1.Caption := '';

Label2.Caption := '';

RadioButton2.Checked := True;

RadioButton2Click(Sender);

LabeledEdit1Change(Sender);

end;


end.


hello_html_m5286a50e.png
















hello_html_45c7732d.png

hello_html_m5a630f1f.png





























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 Lab7

{

public partial class FormSolution : Form

{

private Font drawFontArial10 = new Font("Arial", 10);

private SolidBrush drawBrushLightSkyBlue = new SolidBrush(Color.LightSkyBlue);

private SolidBrush drawBrushBlack = new SolidBrush(Color.Black);


private bool chooseStation1Toggle = false;

private bool chooseStation2Toggle = false;

private bool editRoadsToggle = false;


private PointF station1Point = new PointF(-10, -10);

private PointF station2Point = new PointF(-10, -10);

private Point minRoadToStation1Point;

private Point minRoadToStation2Point;


private Pen drawPenBlue3 = new Pen(Color.Blue, 3);

private Pen drawPenCrimson4 = new Pen(Color.Crimson, 4);


private Point midPoint;


private Point downPoint;

private Point upPoint;


private Point mirrorUpPoint;

private Point mirrorDownPoint;


private int currentAngle = 0;

private int totalDistance=0;


public FormSolution()

{

InitializeComponent();

midPoint.X = 200;

midPoint.Y = 200;

downPoint.X = 0;

downPoint.Y = 0;

upPoint.X = 0;

upPoint.Y = 0;

}


private void btn_exit_Click(object sender, EventArgs e)

{

this.Hide();

}


private void pb_MouseMove(object sender, MouseEventArgs e)

{

reDrawCanvas(new Point(e.X, e.Y));

}


private void reDrawCanvas(Point e)

{

pb.Refresh();

drawStations();

drawAngle();

drawRoads();

drawCurrentCoord(e.X, e.Y);


if ((station1Point.X >= 0) && (station2Point.X >= 0) && !((downPoint.X == 0) && (downPoint.Y == 0) && (upPoint.X == 0) && (upPoint.Y == 0)))

{

calculateAndDrawResult(station1Point, 1);

calculateAndDrawResult(station2Point, 2);

}


if (chooseStation1Toggle == true)

{

drawCircle(e.X, e.Y);

}

else if (chooseStation2Toggle == true)

{

drawCircle(e.X, e.Y);

}


if (editRoadsToggle == true)

{

drawEditRoad(e.X, e.Y);

}

}


private void drawRoads()

{

Graphics g = pb.CreateGraphics();

g.DrawLine(drawPenCrimson4, downPoint, upPoint);

g.DrawLine(drawPenCrimson4, mirrorDownPoint, mirrorUpPoint);


}


private void drawStations()

{

Graphics g = pb.CreateGraphics();


Rectangle drawRectC = new Rectangle(Convert.ToInt32(station1Point.X - 3), Convert.ToInt32(station1Point.Y - 3), 6, 6);

Pen drawPenRed = new Pen(Color.Red, 4);

g.DrawEllipse(drawPenRed, drawRectC);

String StringStation1 = "Пункт 1 (" + Convert.ToInt32(station1Point.X) + ":" + Convert.ToInt32(station1Point.Y) + ")";

g.DrawString(StringStation1, new Font("Arial", 7), drawBrushBlack, station1Point);

drawRectC = new Rectangle(Convert.ToInt32(station2Point.X - 3), Convert.ToInt32(station2Point.Y - 3), 6, 6);

Pen drawPenYellow = new Pen(Color.Yellow, 4);

g.DrawEllipse(drawPenYellow, drawRectC);

String StringStation2 = "Пункт 2 (" + Convert.ToInt32(station2Point.X) + ":" + Convert.ToInt32(station2Point.Y) + ")";

g.DrawString(StringStation2, new Font("Arial", 7), drawBrushBlack, station2Point);


}


private void drawCircle(int X, int Y)

{

Graphics g = pb.CreateGraphics();

Rectangle drawRect = new Rectangle(X - 3, Y - 3, 6, 6);

PointF drawPoint = new PointF(X - 70, Y - 10);

g.DrawEllipse(drawPenBlue3, drawRect);

String drawString;

if (chooseStation1Toggle == true)

{

drawString = "Пункт 1";

}

else

{

drawString = "Пункт 2";

}


g.DrawString(drawString, drawFontArial10, drawBrushBlack, drawPoint);

}

private void drawEditRoad(int X, int Y)

{

Graphics g = pb.CreateGraphics();


Point currentMousePoint = new Point(X, Y);


int upY = 0;

int downY = 400;

int upX=0;

int downX=0;

if (currentAngle <= 90)

{

upY = 0;

upX = lineDrawFormulaX(currentMousePoint, midPoint, upY);


downY = 400;

downX = lineDrawFormulaX(currentMousePoint, midPoint, downY);

}

else if (currentAngle > 90)

{

upX = 0;

upY = lineDrawFormulaY(currentMousePoint, midPoint, upX);


downX = 400;

downY = lineDrawFormulaY(currentMousePoint, midPoint, downX);

}


upPoint = new Point(upX, upY);

downPoint = new Point(downX, downY);


mirrorUpPoint = new Point(400 - upX, upY);

mirrorDownPoint = new Point(400 - downX, downY);

g.DrawLine(drawPenCrimson4, downPoint, upPoint);

g.DrawLine(drawPenCrimson4, mirrorDownPoint, mirrorUpPoint);

}


private void drawAngle()

{

if (!((upPoint.X==0) && (upPoint.Y==0)) || !((downPoint.X==0) && (downPoint.Y==0)))

{

Graphics g = pb.CreateGraphics();


PointF vector1 = new PointF(upPoint.X - downPoint.X, upPoint.Y - downPoint.Y);

PointF vector2 = new PointF(mirrorUpPoint.X - mirrorDownPoint.X, mirrorUpPoint.Y - mirrorDownPoint.Y);


double angle = Math.Acos((vector1.X * vector2.X + vector1.Y * vector2.Y) / (Math.Sqrt(Math.Pow(vector1.X, 2) + Math.Pow(vector1.Y, 2)) * Math.Sqrt(Math.Pow(vector2.X, 2) + Math.Pow(vector2.Y, 2))));

angle = angle * 180 / Math.PI;


currentAngle = Convert.ToInt32(angle);


Rectangle drawRect = new Rectangle(new Point(175,175),new Size(50,50));


int firstCircleAngle =((360 - Convert.ToInt32(angle) * 2) / 4);


g.DrawPie(drawPenBlue3, drawRect, firstCircleAngle, Convert.ToInt32(angle));

g.DrawString(Convert.ToInt32(angle) + "°", drawFontArial10, drawBrushBlack, new PointF(midPoint.X - 15, midPoint.Y + 25));

}



}


private int lineDrawFormulaX(Point po1, Point po2, int currentY)

{

try

{

return (currentY * po1.X - currentY * po2.X - 200 * po1.X + po2.X * po2.Y + po1.Y * po2.X - po2.X * po2.Y) / (po1.Y - po2.Y);

}

catch

{

return (currentY * po1.X - currentY * po2.X - 200 * po1.X + po2.X * po2.Y + po1.Y * po2.X - po2.X * po2.Y) / (po1.Y - 1 - po2.Y);

}

}


private int lineDrawFormulaY(Point po1, Point po2, int currentX)

{

try

{

return (po2.Y * currentX - po2.Y * po1.X - po1.Y * currentX + po1.Y * po1.X + po1.Y * po2.X - po1.Y * po1.X) / (po2.X - po1.X);

}

catch

{

return (po2.Y * currentX - po2.Y * po1.X - po1.Y * currentX + po1.Y * po1.X + po1.Y * po2.X - po1.Y * po1.X) / (po2.X - 1 - po1.X);

}

}



private void drawCurrentCoord(int X, int Y)

{

Graphics g = pb.CreateGraphics();

PointF drawPoint = new PointF(X, Y - 20);

String drawString = "(" + Convert.ToString(X) + " : " + Convert.ToString(Y) + ")";

g.DrawString(drawString, drawFontArial10, drawBrushLightSkyBlue, drawPoint);

}


private void btn_choose_station_1_Click(object sender, EventArgs e)

{

chooseStation1Toggle = true;

chooseStation2Toggle = false;

editRoadsToggle = false;

btn_choose_station_1.Enabled = false;


}


private void pb_MouseClick(object sender, MouseEventArgs e)

{

if (chooseStation1Toggle == true)

{

chooseStation1Toggle = false;

chooseStation2Toggle = false;

editRoadsToggle = false;

btn_choose_station_1.Enabled = true;

btn_choose_station_2.Enabled = true;

station1Point.X = e.X;

station1Point.Y = e.Y;

}

else if (chooseStation2Toggle == true)

{

chooseStation1Toggle = false;

chooseStation2Toggle = false;

editRoadsToggle = false;

btn_choose_station_1.Enabled = true;

btn_choose_station_2.Enabled = true;

station2Point.X = e.X;

station2Point.Y = e.Y;

}

else if (editRoadsToggle == true)

{

chooseStation1Toggle = false;

chooseStation2Toggle = false;

editRoadsToggle = false;

btn_choose_station_1.Enabled = true;

btn_choose_station_2.Enabled = true;

}


}


private void btn_choose_station_2_Click(object sender, EventArgs e)

{

chooseStation2Toggle = true;

chooseStation1Toggle = false;

editRoadsToggle = false;

btn_choose_station_1.Enabled = false;

}


private void FormSolution_Load(object sender, EventArgs e)

{

}


private void pb_Click(object sender, EventArgs e)

{


}


private void button1_Click(object sender, EventArgs e)

{

editRoadsToggle = true;

chooseStation1Toggle = false;

chooseStation2Toggle = false;

}


private void btn_calculate_Click(object sender, EventArgs e)

{

calculateAndDrawResult(station1Point, 1);

calculateAndDrawResult(station2Point, 2);


}


private void calculateAndDrawResult(PointF Station, int Key)

{


int Y, X = 0, D, i = 0;

int[,] distanceAndXYArray = new int[3, 802];


if (((downPoint.X == 400) && (downPoint.Y == 200)) || ((downPoint.X == 0) && (downPoint.Y == 200)))

{

for (X = 0; X < 401; X++)

{

Y = lineDrawFormulaY(downPoint, midPoint, X);

D = getLengthOfTwoPoints(new Point(Convert.ToInt32(Station.X), Convert.ToInt32(Station.Y)), new Point(X, Y));

distanceAndXYArray[0, i] = X;

distanceAndXYArray[1, i] = Y;

distanceAndXYArray[2, i] = D;

i++;

}


for (X = 0; X < 401; X++)

{


Y = lineDrawFormulaY(mirrorDownPoint, midPoint, X);

D = getLengthOfTwoPoints(new Point(Convert.ToInt32(Station.X), Convert.ToInt32(Station.Y)), new Point(X, Y));

distanceAndXYArray[0, i] = X;

distanceAndXYArray[1, i] = Y;

distanceAndXYArray[2, i] = D;

i++;

}

}

else

{

for (Y = 0; Y < 401; Y++)

{

X = lineDrawFormulaX(downPoint, midPoint, Y);

D = getLengthOfTwoPoints(new Point(Convert.ToInt32(Station.X), Convert.ToInt32(Station.Y)), new Point(X, Y));

distanceAndXYArray[0, i] = X;

distanceAndXYArray[1, i] = Y;

distanceAndXYArray[2, i] = D;

i++;

}


for (Y = 0; Y < 401; Y++)

{


X = lineDrawFormulaX(mirrorDownPoint, midPoint, Y);

D = getLengthOfTwoPoints(new Point(Convert.ToInt32(Station.X), Convert.ToInt32(Station.Y)), new Point(X, Y));

distanceAndXYArray[0, i] = X;

distanceAndXYArray[1, i] = Y;

distanceAndXYArray[2, i] = D;

i++;

}

}



int minElement = distanceAndXYArray[2, 0];

int minI = 0;

for (i = 0; i < 800; i++)

{

if (minElement > distanceAndXYArray[2, i])

{

minElement = distanceAndXYArray[2, i];

minI = i;

}

}


Graphics g = pb.CreateGraphics();

Rectangle drawRect;

if (Key == 1)

{

minRoadToStation1Point.X = distanceAndXYArray[0, minI];

minRoadToStation1Point.Y = distanceAndXYArray[1, minI];


drawRect = new Rectangle(minRoadToStation1Point.X - 3, minRoadToStation1Point.Y - 3, 6, 6);

g.DrawEllipse(new Pen(Color.Black, 2), drawRect);

g.DrawLine(new Pen(Color.Black, 2), Station, minRoadToStation1Point);

}

else if (Key == 2)

{

minRoadToStation2Point.X = distanceAndXYArray[0, minI];

minRoadToStation2Point.Y = distanceAndXYArray[1, minI];

drawRect = new Rectangle(minRoadToStation2Point.X - 3, minRoadToStation2Point.Y - 3, 6, 6);

g.DrawEllipse(new Pen(Color.Black,2), drawRect);

g.DrawLine(new Pen(Color.Black, 2), Station, minRoadToStation2Point);


}


totalDistance = 0;

totalDistance += getLengthOfTwoPoints(new Point(Convert.ToInt32(station1Point.X), Convert.ToInt32(station1Point.Y)), minRoadToStation1Point);

totalDistance += getLengthOfTwoPoints(new Point(Convert.ToInt32(station2Point.X), Convert.ToInt32(station2Point.Y)), minRoadToStation2Point);

totalDistance += getLengthOfTwoPoints(minRoadToStation1Point, midPoint);

totalDistance += getLengthOfTwoPoints(minRoadToStation2Point, midPoint);

g.DrawString("Расстояние пути: " + Convert.ToString(totalDistance), drawFontArial10, drawBrushBlack, new PointF(50, 20));


}


private int getLengthOfTwoPoints(Point po1, Point po2)

{

return Convert.ToInt32(Math.Sqrt(Math.Pow((po2.X-po1.X),2)+Math.Pow((po2.Y-po1.Y),2)));


}


private void pb_Paint(object sender, PaintEventArgs e)

{


}

}

}





























2.10 Лабораторная работа 10

Тема. Моделирование объектов методом

пространства состояния, динамика которого

описывается дифференциальным уравнением


Задачи исследования:

Исследовать динамику объекта методом пространства состояния с использованием информационных технологий MatLab методом составления программы: решения уравнений пространства состояний.

Исходные данные

Начальные и конечные координаты объекта x(0), x`(0), xk, такт исследования T, управляющее воздействие u(0). Получить зависимости x(T), x`(T) с использованием

программ MatLab , Scilab

Математическая модель объекта.

Исходное дифференциальное уравнение имеет вид

hello_html_m2e2f8e7a.gif


Представим исходное уравнение в нормальной форме Коши, вводя обозначения

hello_html_m5f9d7872.gif


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

hello_html_ed9157c.gif



С целью определения матриц уравнения параметров состояния преобразуем систему:



hello_html_m7d01ad13.gif



матрицы коэффициентов соответственно равныhello_html_m7d01ad13.gifhello_html_7017e8e6.gifhello_html_m705c5775.gifhello_html_m37f07e0b.gifhello_html_18ed2da.gif


Уравнение параметров состояния будет иметь вид


hello_html_m157ec255.gif




Решение данных примеров в системе Matlab

Для этого создадим 3 m-файла (для каждого примера отдельно). Далее эти m-файлы запустим в командном окне(command Window).

Описание m-файлов

Пример 1

figure(1); //вызов графического окна;

clf; //очистка фигуры;

title(“solve equation d(dx)+k1*dx=k2*u”); //установка титульной надписи

syms k1 k2 p t T i Y01 Y02; //задание переменных;

//ввод данных;

k2=0.01;

k1=0.5;

//ввод исходных матриц;

B=[0 1;0 -k1];

U=[0.1];

E=eye(size(B)); //создание единичной матрицы

A=[0;k2];

D=E*p;

W=D-B;

W1=inv(W); //нахождение обратной матрицы

W2=ilaplace(W1,p,t); //обратное преобразование Лапласа

H=(int(W2,t,0,T))*A; //матрица управления

Y0=[Y01;Y02]; //матрица начальных условий

Y01=0;Y02=1;

W21=ilaplace(W1,p,T); //матрица состояния

T=0.1; // такт работы

for i=0:10,

Y=W21*Y0+H*U; //расчетная формула

Y1=Y01+(-2*exp(-1/2*T)+2)*Y02+1/500*T+1/250*exp(-1/2*T)-1/250;

Y2=exp(-1/2*T)*Y02-1/500*exp(-1/2*T)+1/500;

Y01=Y1; Y02=Y2;

hold on;.//обеспечение продолжение вывода графиков в текущее окно

plot(i,Y01,’o’,i,Y02,’or’);//построение графиков

legend(“x”,’dx’); //добавление к текущему графику легенды

xlabel(“i”); //установка надписи на оси OX

ylabel(“x,dx”); //становка надписи на оси OYgrid on; //сетка

end

Y2=exp(-1/2*T)*Y02-1/500*exp(-1/2*T)+1/500;

Y01=Y1; Y02=Y2;

hold on;.//обеспечение продолжение вывода графиков в текущее окно

plot(i,Y01,’o’,i,Y02,’or’);//построение графиков

legend(“x”,’dx’); //добавление к текущему графику легенды

xlabel(“i”); //установка надписи на оси OX

ylabel(“x,dx”); //становка надписи на оси OYgrid on; //сетка

end



нач. Условия


k1=0.5, k2=0.01; i=0..10; y10=0; y20=1; n=0.1; u=0.1





hello_html_4552311.png




Y2=exp(-1/2*T)*Y02-1/500*exp(-1/2*T)+1/500;

Y01=Y1; Y02=Y2;

hold on;.//обеспечение продолжение вывода графиков в текущее окно

plot(i,Y01,’o’,i,Y02,’or’);//построение графиков

legend(“x”,’dx’); //добавление к текущему графику легенды

xlabel(“i”); //установка надписи на оси OX

ylabel(“x,dx”); //становка надписи на оси OYgrid on; //сетка

end

hello_html_ff60b12.gif



Пример. Для дифференциального уравнения

hello_html_m59ca5fcc.gif


Сhello_html_m734071c7.gifhello_html_28146482.gifhello_html_26048a56.gif целью определения матриц уравнения параметров состояния преобразуем систему:

hello_html_m15c900f0.gifhello_html_m7d28e756.gifhello_html_m7c8fd428.gifhello_html_m7c8fd428.gifhello_html_m657ecfed.gifhello_html_68297d8a.gifhello_html_1a2a11b8.gifhello_html_332dbd81.gif






Пример 2

figure(2); //вызов графического окна;

clf; //очиска фигуры;

title(“solve equation d(dx)+2*dx-3*x=u”); //установка титульной надписи

syms p t T i Y10 Y20; //задание переменных;

//ввод данных;

B=[0 1;3 -2];

E=eye(size(B)); //создание единичной матрицы

U=[0.1];

A=[0;1];

D=E*p;

W=D-B;

W1=inv(W); //нахождение обратной марицы

W2=ilaplace(W1,p,t); //обратное преобразование Лапласа

H=(int(W2,t,0,T))*A; //марица управления

Y10=0;

Y20=1;

T=0.1; // такт работы

for i=0:10,

Y=W21*Y0+H*U;

Y1=(1/4*exp(-3*T)+3/4*exp(T))*Y10+1/16*16^(1/2)*(exp(T*(-

1+1/2*16^(1/2)))-exp(T*(-1-1/2*16^(1/2))))*Y20+1/40*(exp(-

T)*exp(1/2*16^(1/2)*T)*16^(1/2)+8*exp(-T)*exp(1/2*16^(1/2)*T)-exp(-

T)*exp(-1/2*16^(1/2)*T)*16^(1/2)+8*exp(-T)*exp(-1/2*16^(1/2)*T)-

16)/(-2+16^(1/2))/(2+16^(1/2));

Y2=3/16*16^(1/2)*(exp(T*(-1+1/2*16^(1/2)))-exp(T*(-1-

1/2*16^(1/2))))*Y10+(3/4*exp(-3*T)+1/4*exp(T))*Y20-1/40*exp(- 3*T)+1/40*exp(T);

Y10=Y1;

Y20=Y2;

hold on; //обеспечение продолжение вывода графиков в текущее окно

plot(i,Y1,’o’,i,Y2,’or’); //построение графиков

legend(“x”,’dx’); //добавление к текущему графику легенды

xlabel(“i”); //установка надписи на оси OX

ylabel(“x,dx”); //становка надписи на оси OY

grid on; end //сетка

hello_html_mbb1ea7c.gif





















3. Индивидуальные задания

по моделированию

Целью работ по компьютерному моделированию являются:

- выработка практических навыков в выполнении работ по математическому моделированию;

- освоение элементов самостоятельной научно- исследовательской работы;

- закрепление навыков программирования, полученных на занятиях.

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

* постановка задачи ( начальные условия);

* математическая модель;

* метод исследования модели;

* алгоритм моделирования задачи ( блок-схема);

* программа на языке программирования, заданным преподавателем;

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

( табличный, графический, динамический);

* содержательный анализ результатов моделирования, выводы. 3.1 Движение тел в среде с учетом трения

1. Парашютист прыгает с некоторой высоты и летит, не открывая парашюта; на какой высоте и через какое время ему следует открыть парашют, чтобы к моменту приземления иметь безопасную скорость (не большую 10 м/с)?

2 Исследовать, как связана высота прыжка с площадью поперечного сечения пара­шюта, чтобы скорость приземления была безопасной.(10 м/с)

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

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

5. Промоделировать движение исследовательского зонда, «выстреленного» вертикально вверх с уровня земли. В верхней точке траектории над зондом раскрывается парашют, и он плавно спускается в точку старта. Исследовать влияние ветра.

6. Промоделировать движение исследовательского зонда, «выстреленного» верти­кально вверх с летящего над землей самолета. В верхней точке траектории над зондом раскрывается парашют, и он плавно спускается на землю. Исследовать влияние ветра.

7. Глубинная бомба, установленная на взрыв через заданное время, сбрасывается со стоящего неподвижно противолодочного корабля. Исследовать связь между глубиной, на которой произойдет взрыв, и формой корпуса (сферической, полусфе­рической, каплевидной и т.д.).

­8. Глубинная бомба, установленная на взрыв на заданной глубине, сбрасывается со стоящего неподвижно противолодочного корабля. Исследовать связь между вре­менем достижения заданной глубины и формой корпуса (сферической, полусфе­рической, каплевидной и т.д.).

9. Провести моделирование взлета ракеты при значениях параметров m0 ,m(кон) ,Fтяги ,a. При каких параметрах ракета может достигнуть значения первой космической скорости.

10. Провести исследование соотношения входных параметров m0 и Fтяги, при которых ракета достигнет первой космической скорости (и в соответствующий момент исчерпает горючее). Остальные входные параметры фиксировать произвольно. По­строить соответствующую фазовую диаграмму в переменных (mo и Fтяги).

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

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

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

14. Глубинная бомба-торпеда, снабженная разгонным двигателем, установленная на взрыв через заданное время, сбрасывается со стоящего неподвижно противоло­дочного корабля. Исследовать связь между глубиной, на которой произойдет взрыв, и формой корпуса (сферической, полусферической, каплевидной и т.д.).

15 Глубинная бомба-торпеда, снабженная разгонным двигателем, установленная на взрыв на заданной глубине, сбрасывается со стоящего неподвижно противоло­дочного корабля. Исследовать связь между временем достижения заданной глубины, и формой корпуса (сферической, полусферической, каплевидной).

16. Глубинная бомба-торпеда, снабженная разгонным двигателем, нацеливается с подводной лодки на стоящий неподвижно противоло­дочный корабль. Исследовать связь между временем достижения заданной глубины и формой корпуса (сферической, полусферической, каплевидной и т.д.).

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

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

19. Поставить и решить задачу о подводной охоте при дополнительном условии: акула движется.

20. Промоделировать движение исследовательского зонда, «выстреленного» под углом к горизонту. В верхней точке траектории над зондом раскрывается тормозной парашют, затем зонд плавно движется до земли.

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

22. Глубинная бомба-торпеда, снабженная разгонным двигателем, установлен» на взрыв на заданной глубине, сбрасывается с движущегося противолодочни корабля. Исследовать связь между временем достижения заданной глубины, прои-денным расстоянием по горизонтали и формой корпуса (сферической, полусфе­рической, каплевидной и т.д.).

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

24 Торпеда, снабженная разгонным двигателем, нацеливается с подводной лодки на стоящий вертикально над ней надводный корабль. Исследовать связь между временем поражения цели и формой корпуса (сферической, полусферической, каплевидной и т.д.).

25 Построить траектории и найти временные зависимости горизонтальной и вер­тикальной составляющих скорости и перемещения для тела массой «м» кг, брошен­ного под углом «а» к горизонту с начальной скоростью «v» м/с:1) в воздухе; 2) в воде.

26. Смоделировать полет камня без учета силы трения. Камень массой « m» брошен под углом «а» к горизонту со скоростью «v». Представить траекторию полета без учета сопротивления воздуха. Исследовать как ме­няются максимальная высота и дальность полета камня при изменении угла . Показать графически.

27. Смоделировать полет бумажки без учета силы трения. Исследовать влияние массы на дальность и высоту полета? Результаты представить графически.

28.Смоделировать полет камня с учетом сопротивления воздуха и исследовать влияние коэффициентов на максимальные дальность и высоту полета. Камень массой «m» брошен под углом «a» к горизонту со скоростью «v» м/с. Cравнить траектории полета без учета сопротивления воздуха.

29. Дальность полета. Как меняются максимальная дальность полета и время полета комка бумаги массой 20 г в зави­симости от угла бросания?

30. Оптимальный угол бросания бумажки. Найдите оптималь­ный угол бросания комка бумаги для получения максимальной дальности полета.

31. Оптимальный угол бросания камня. Под каким углом к го­ризонту следует бросить камень массой 200 г со скоростью 20 м/с, чтобы дальность полета была наибольшей с учетом силы сопротивле­ния воздуха Сравните со случаем, когда сопротивление воздуха не учитывается.

32. Смоделировать полет бумажки без учета силы трения. Исследовать влияние массы на дальность и высоту полета? Результаты представить графически.

33. Смоделировать полет камня с учетом сопротивления воздуха и исследовать влияние коэффициентов на максимальные дальность и высоту полета. Камень массой «m» брошен под углом «a» к горизонту со скоростью «v» м/с. Cравнить траектории полета без учета сопротивления воздуха.

34. Исследовать полет комка бумаги. Как меняются максимальная дальность полета и время полета комка бумаги массой «m» «г» в зави­симости от угла бросания? Начальная скорость комка бумаги рав­на «v» «м/с», сила трения «F».

35. Определить оптимальный угол бросания бумажки массой «m» «г», начальная скорость рав­на «v» «м/с», сила трения «F». 36. Исследовать приближение космического аппарата к Луне. Космический объект массой «m» кг подлетает к Луне. Когда расстояние «а» становится равным «s» км и прицельное расстояние «р» км и скорость объекта «v»м/с.

37. Рассчитайте траекторию полета вблизи Лу­ны. Масса Луны равна 7.3*10^22 кг, радиус Луны равен 1.7*10^3 км. Напишите уравнения движения по оси х и по оси у с учетом зависимости силы притяжения от расстояния.

38. Cмоделировать траектории подлета космического аппарата вблизи Луны. Пусть скорость объекта равна v перпендикулярна направлению на центр Луны. Прицельное расстоя­ние p=5*10^3 км. Как будет двигаться космический аппарат при скоростях 500 м/с, 1000 m/c, 1500 m/c. Радиус Луны ранен 1.7*10^3 км

39. Cмоделировать различные космические траектории. За­дать различные значения для величии р, а, v и исследовать вопрос о движении аппарата вблизи Луны (или Земли) более подробно.

9. Cмоделировать посадку спутника в атмосфере. Спутник Земли массой «м» «т» вошел в атмосферу на высоте S км со скоростью V км/с, параллельной поверхности Земли. Атмосфера тормозит полет силой, равной Av, где для простоты расчетов выражение А/т не зависит от высоты и равно 10^(-4) с-1 .Начертите траекторию посадки. Радиус Земли равен 6370 км, масса Земли равна 5,96*10^24 кг.

10. Исследовать движение болида. Каменная глыба (болид) массой «м» «т» приближает­ся из космоса к планете, у которой радиус R= 1,74*106 м. ,масса т == 7,3*1022 кг и толщина атмосферы h=1,26*106 м. Пусть трение в атмосфере характеризуется силой Av, причем A/m=5*10^(-4) c-1. (коэффициент А не зависит от высоты).

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

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

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

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

15. Как будет выглядеть полет искусственного спутника Земли, если учесть возму­щающее действие Луны?

16. Разработать и реализовать модель движения искусственного спутника Земли при учете воздействия на него малой постоянной силы, обусловленной «солнеч­ным ветром». Считать, что плоскость орбиты движения спутника изначально пер­пендикулярна к «солнечному ветру».

17. Считая, что движение Луны вокруг Земли происходит практически по круго­вой орбите, проанализировать воздействие на эту орбиту со стороны Солнца для малого участка движения, на котором плоскость орбиты перпендикулярна к оси «Солнце—Земля».

18. Проанализировать особенности движения искусственного спутника Земли, дви­жущегося практически по круговой орбите на высоте порядка 300 км, связанные с малым сопротивлением атмосферы.

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

Моделирование биологических систем

Модель однородной популяции

X(i+1) = x(i)+ax(i)-bx(i)2; x0 = c, i = 0.1...n, (3.82)

где n предельное время моделирования.

Усложнение модели может быть выполнено за счет учета переменности коэффициентов - a (скорости роста) и b (скорости гибели) - в зависимости от времени, а также с учетом половых, возрастных различий индивидумов.

Модель межвидовой конкуренции


dN1/dt=r1*N1(K1-N1-a12*N2)/k1

dN2/dt=r2*N2(K2-N2-a21*N1)/k2

где k1,k2 плотности насыщения

r1, r2 врожденные скорости роста

a12, a21- коэффициенты конкуренции



Эпидемия болезней

В изолированном поселке с населением m человек возникла эпидемия болезни, распространение которой описывается соотношениями:

xi+1 = xi – b*xi*yi;

yi+1 = yi - cyi + bxi*yi;

zi+1 = zi + cyi;

x0=a0, y0=b0, z0=c0,

где xi, yi, zi - число здоровых, больных (инфицированных) и невосприимчивых (переболевших) в момент времени i=0.1...n;

b - частота контактов больных и здоровых;

c - величина, обратная среднему времени выздоровления и зависящая от эффективности лекарств 0

Более строго соотношение может быть получено из системы дифференциальных уравнений:


Модель “хищник - жертва”

Имеются популяции двух видов, которые представляются отношениями

xi+1 = xi + a1xi - b1xi2 - g1xiyi , x0 = c1,

yi+1 = yi - a2yi + b2yi2 - g2xiyi , y0 = c2,

где xi - численность (плотность) жертв,

yi - численность хищников,

g1 - коэффициент защиты жертв,

g2 - коэффициент прожорливости хищников.


Рост опухоли

Раковая опухоль обычно увеличивается экспоненциально в соответствии с дифференциальным уравнением:

hello_html_m17d3c5fa.gif

, где v - размер опухоли,

с,a,b - константы.

Определить, при каких значениях параметров С существует предельный размер опухоли. Выяснить, при каких значениях С рост опухоли не превосходит некоторой конечной величины.

hello_html_0.gifAutoShape 16Frame1AutoShape 9


1. Изучить характер эволюции популяции, при зна­чениях параметров a, b, n в зависимости от значения параметра b в диапазоне 0,1 < b < 10,

Есть ли качественные различия в характере эволюции в зависимости от значе­ния b?

2. Изучить характер эволюции популяции при зна­чениях параметров b, a, n в зависимости от значения параметра a в диапазоне 1 < a < 10,

Есть ли качественные различия в характере эволюции в зависимости от значе­ния a?

3. Изучить характер эволюции популяции а, b,v в зависимости от значения параметра V в диапазоне 1 <v<10.

Есть ли качественные различия в характере эволюции в зависимости от значе­ния v?

4. Реализовать модель при различных наборах значений параметров: n, a, b и изучить влияние параметров на характер эволюции.

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

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

7. Реализовать моделирование межвидовой конкуренции по формулам при значениях параметров r1 = 2,r2= 2,k1 = 200, k2 = 200, a12 = 0,5,a21= 0,5. Проана­лизировать зависимость судьбы популяций от соотношения значений их началь­ной численности n1,n2.

8. Реализовать моделирование межвидовой конкуренции по формулам при значениях параметров r1 = 2,r2 = 2, k1 = 200, k2= 200, a12= 100, a21 = 100. Проанализировать зависимость судьбы популяций от соотношения значений ко­эффициентов конкуренции a12, a21,.

9. Построить в фазовой плоскости (n1n2) границы зон, разделяющих какие-либо два режима эволюции конкурирующих популяций (в соответствии с моделью ). Остальные параметры модели выбрать произвольно. Учесть при этом, что режим устойчивого сосуществования популяций может в принципе реализовать только при ф12 < 1.

10. Провести моделирование динамики численности популяций в системе «хищник-жертва» при значениях параметров r = 5, a= = 0, q = 2 c= 0,6 Проанализировать зависимость исхода эволюции от соотношения значений пара­метров n, c0.

11. Провести моделирование динамики численности популяций в системе «хищ­ник—жертва» при значениях параметров r = 5, а = 0,1, q = 2, f = 100, Со = 6. Проанализировать зависимость результатов моделирования от значения параметра a в диапазоне 0,1 : 2.

12. Провести моделирование динамики численности популяций в системе «хищ­ник—жертва» при значениях параметров r = 5, а = 0,1, f= 2, N = 100,

13. Провести моделирование динамики численности популяций в системе «хищ­ник—жертва» при значениях параметров а = 0,! f= 2, q = 2, N= 100, Со=б. Проанализировать зависимость результатов моделирования от значе­ния параметра r в диапазоне 0,1 < r < 2.

14. Провести моделирование динамики численности популяций в системе «хищ­ник—жертва» при значениях параметров a, q, No ,Co. Проанализировать зависимость результатов моделирования от значе­ния параметра «а» в диапазоне 0,1 < а < 2.

15. Модель «хищ­ник—жертва» предсказывает сопряженные колебания численности жертв и хищников. Исследовать зависимость запаздывания амплитуд колебаний числен­ности хищников от амплитуд колебаний численности жертв в зависимости от значений параметра а. Значения остальных параметров фиксировать по усмотре­нию.

16. Модель «хищ­ник—жертва» предсказывает сопряженные колебания численности жертв и хищников. Исследовать зависимость запаздывания амплитуд колебаний числен­ности хищников от амплитуд колебаний численности жертв в зависимости от значений параметра q. Значения остальных параметров фиксировать по усмотре­нию.

17.Модель предсказывает сопряженные колебания численности жертв и хищников. Исследовать зависимость запаздывания амплитуд колебаний численно­сти хищников от амплитуд колебаний численности жертв и зависимости от значе­ний параметра a. Значения остальных параметров фиксировать по усмотрению.

18. Модель предсказывает сопряженные колебания численности жертв и хищников. Исследовать зависимость запаздывания амплитуд колебаний численно­сти хищников от амплитуд колебаний численности жертв в зависимости от значе­ний параметра f, Значения остальных параметров фиксировать по усмотрению.

19.Модель предсказывает сопряженные колебания численности жертв и хищников. Исследовать зависимость запаздывания амплитуд колебаний численно­сти хищников от амплитуд колебаний численности жертв в зависимости от соот­ношения значений начальных численностей популяции N0, С0. Значения осталь­ных параметров фиксировать по усмотрению.

2hello_html_28311771.gif3.Промоделировать процесс распространения инфекции стригущего лишая по участку кожи размером n*n (п — нечетное) клеток.

1.Внутривидовая конкуренция в популяции с


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

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

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

Моделирование случайных процессов

1. Провести моделирование очереди в магазине с одним продавцом при равнове­роятных законах распределения описанных выше случайных величин: прихода по­купателей и длительности обслуживания (при некотором фиксированном наборе параметров). Получить устойчивые характеристики: средние значения ожидания в очереди покупателем и простой продавца в ожидании прихода покупателей. Оце­нить их достоверность. Оценить характер функции распределения величин g и h.

2. Провести то же моделирование при пуассоновских законах распределения ве­роятностей входных событий: прихода покупателей и длительности обслуживания (при некотором фиксированном наборе параметров).

3.Провести то же моделирование при нормальном законе распределения вероят­ностей входных событий: прихода покупателей и длительности обслуживания (при некотором фиксированном наборе параметров).

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

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

6.Смоделировать ситуацию, описанную в предыдущем варианте, но считать, что, если в момент попытки сделать заказ обе телефонистки занять;, формируется очередь.

7. Пусть на телефонной станции с одним входом используется обычная система: если абонент занят, то очередь не формируется и надо звонить снова. Смоделиро­вать ситуацию: три абонента пытаются дозвониться до одного и того же владельца номера и в случае успеха разговаривают с ним некоторое (случайное по длитель­ности) время. Какова вероятность того, что некто, пытающийся дозвониться, не сможет сделать это за определенное время Т?

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

9. На травм. пункте работает один врач. Длительность лечения больного и проме­жутки времени между поступлениями больных — случайные величины, распреде­ленные по пуассоновскому закону. По тяжести травм больные делятся на три кате­гории, поступление больного любой категории — случайное событие с равноверо­ятным распределением. Врач вначале занимается больными с максимально тяже­лыми травмами (в порядке их поступления), затем, если таковых нет, — больны­ми с травмами средней тяжести (в порядке их поступления) и лишь затем — боль­ными с легкими травмами. Смоделировать процесс и оценить средние времена ожидания в очереди больных каждой из категорий,

10.Смоделировать ситуацию, описанную в предыдущем варианте, при условии, что в травм. пункте работают два врача, а больные делятся не на три, а на две категории.

11.Одна ткачиха обслуживает группу станков, осуществляя по мере необходимос­ти краткосрочное вмешательство, длительность которого — случайная величина.Какова вероятность простоя сразу двух станков? Как велико среднее время про­стоя одного станка?

12.Смоделировать ситуацию, описанную в предыдущем варианте, если группу стан­ков совместно обслуживают две ткачихи.

13.В городском автохозяйстве две ремонтные зоны. Одна — обслуживает ремонты краткой и средней продолжительности, другая — средней и долгой (т.е. средне­срочный ремонт может осуществлять каждая из зон). По мере поломок в автохо­зяйство доставляют транспорт; промежуток времени между доставками — случай­ная пуассоновская величина. Продолжительность ремонта — случайная величина с нормальным законом распределения. Смоделировать описанную систему. Каковы средние времена ожидания в очереди транспорта, требующего соответственно крат­косрочного, среднесрочного и длительного ремонта?

14. Реализовать имитационную модель статистического моделирования для реше­ния задачи Бюффона (XVIII в.). Автор аналитически нашел, что если на поле, разграфленное параллельными прямыми, расстояние между которыми L, броса­ется наугад игла длиной l, то вероятность того, что игла пересечет хотя бы одну прямую, определяется формулой р = 2*l/(pi*L) .

Эта задача дала способ имитационному определению числа pi. Действительно, если L =2*l, то p = 1/pi. В ходе моделирования выполнить этот расчет.

15. Разработать модель случайного одномерного блуждания (модель «пьяницы»). Блуждание задается по правилу: если случайное число из отрезка [0,1] меньше 0,5, то делается шаг вправо на расстояние h, в противном случае ~ влево. Распределе­ние случайных чисел принять равновероятным. Решить задачу: какова вероятность при таком блуждании удалиться от началь­ной точки на п шагов?

16.В условиях задачи из предыдущего варианта получить ответ на вопрос: какова вероятность «пьяницей вернуться через п шагов в начальную точку?

17.Точка хаотически блуждает на плоскости по узлам квадратной сетки с возмож­ностью делать с равной вероятностью шаги влево-вправо-вверх-вниз на фиксиро­ванный (за один ход) шаг. Движение происходит в замкнутом прямоугольном объе­ме, и при соприкосновении со стенкой происходит зеркальное отражение от нее.

Ответить в ходе моделирования на вопрос: как связана частота посещения каж­дого узла с расстоянием от него до того узла, из которого начинается движение?

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

19.Смоделировать полет пчелы. На плоскости (поляне) случайным образом растут медоносные растения с заданной концентрацией (на 1 м3). В центре — улей, из кото­рого вылетает пчела. Пчела может долететь от одного растения до любого другого растения, но вероятность выбора монотонно уменьшается с увеличением расстояния между растениями (по некоторому закону). Какова вероятность посещения пчелой конкретного заданного растения за заданное количество элементарных перелетов?

20. Реализовать модель плоского броуновского движения п частиц в прямоугольни­ке. Частицы считать шариками конечного размера. Удары частиц друг о друга и о стенки моделировать как абсолютно упругие. Определить п этой модели зависи­мость давления газа на стенки от числа частиц.

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

22. Реализовать имитационную модель системы «хищник—жертва» по следующей схеме.

«Остров» размером 20^20 заселен дикими кроликами, волками и волчицами, Имеется по несколько представителей каждого вида. Кролики в каждый момент перемещаются с одинаковой вероятностью 1/9 передвигаются в один из восьми сосед­них квадратов (за исключением участков, ограниченных береговой линией) или просто сидят неподвижно. Каждый кролик с вероятностью 0,2 превращается в двух кроликов. Каждая волчица передвигается случайным образом, пока в одном из соседних восьми квадратов не окажется кролик, за которым она охотится. Если волчица и кролик оказываются в одном квадрате, волчица съедает кролика и по­лучает одно очко. В противном случае она теряет 0,1 очка.

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

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

23.Промоделировать процесс распространения инфекции стригущего лишая по участку кожи размером n*n (п — нечетное) клеток.



3.5 Моделирование оптимальных систем


1 Cмоделировать маршрут движения катера Внутри водоема правильной круглой формы радиуса R расположен маленький островок радиуса r. Вычислите и укажите кратчайший прямой марш­рут катера, соединяющий какие-нибудь точки берега и имеющий промежуточный причал у островка.

2 Место для завода

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

.3 Газетный киоск

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

4 Где построить школу?

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

5. С наименьшей суммой расстояний

Три населенных пункта расположены в вершинах ост­роугольного треугольника. Где нужно построить завод, чтобы сумма расстояний от него до всех трех данных пунк­тов была иаименьшей?

6. Проселочная дорога

Через город проходит магистраль, на некотором рассто­янии от которой находится населенный пункт.

7. Направление магистрали

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

8. Наилучшее расположение

Как должна проходить магистраль, чтобы сумма рас­стоянии от нее до трех данных населенных пунктов была наименьшей?

9. Выбор маршрута

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

Каким должен быть кратчайший замкнутый маршрут автобуса, предназначенного для развозки жителей насе­ленного пункта по всем трем заводам?

10. Как проложить дорогу?

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

11. Кратчайший замкнутый маршрут

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

12. Строительство водопровода

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

13. Кратчайшая дорога

Магистраль п канал пересекаются под углом меньше 45°, внутри которого расположен населенный пункт. Как проложить кратчайшую дорогу, проходящую от одного пункта сначала к берегу канала, а затем к магистрали?

14. Мост через какал

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

15. Железнодорожная платформа

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

16. Кратчайший маршрут

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


Задания для метода « Линейное программирование»




hello_html_6dd4bb91.png




Индивидуальные задания

- модель движения материальной точки Аристотеля и Ньютона;

- модель Солнечной системы Птолемея, Коперника, Кеплера;

- простейшую демографическую модель;

- модель многоотраслевой экономики Леонтьева;

- модель процесса распространения эпидемий;

- модель динамики численности биологических популяций;

- модель относительных движений в классической механике;

- модель остывания нагретых тел в атмосфере;

- математическую модель процесса загрязнения воды;

- модель колебательных процессов в физике.

- модель биологической системы «хищник-жертва»;

- модель биологической системы конкурирующих популя­ций;

- модель поведение динамики многочастичной системы;

- модели марковских случайных процессов;

- модель поведение динамической системы, описываемой

уравнени­ями Колмогорова.;

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

- модель исследования систем на основе матричных методов пространства состояния;

- модели экспертных систем;

- модели управления педагогическими системами;

- модели геоинформационных систем .












ЛИТЕРАТУРА


1. Алексеев, О. Г. Управление в системах РАВ / О. Г. Алексеев и др. - М. : МО, 1980. - 600 с.

2. Белошапка, В. Информатика как наука о буквах / В. Белошапка // Информатика и образование. - 1992. - № 1. -С 6-12.

3. Вершинин, О. Е. Компьютер для менеджера / О. Е. Вершинин. - М. : Высшая школа, 1990. - 240 с.

4. Глушков, В. М. Кибернетика / В. М. Глушков. - М. : Наука, 1986.

5. Глушков, В. М. Введение в АСУ / В. М. Глушков. -Киев. : Техника, 1974.

6. Глушков, В. М. Основы безбумажной информатики / В. М. Глушков. - М. : Наука, 1978.

7. Горстко, А. Б. Познакомтесь с математическим моделированием / А. Б. Горстко. - М. : Знание, 1991.

8. Гулд, Х. Компьютерное моделирование в физике /Х. Гулд и др. - М. : Мир, 1990.

9. Ершов, А. П. Об информационной модели машин / А. П. Ершов // Микропроцессоры, средства и системы. - 1985. - № 4. - С 2. 10. Заварыкин, В. М. Численные методы / В. М. Заварыкин. - М. : Просвещение, 1991.

11. Ивахненко, А. Т. Перцептрон - система распознавания образов / А. Т. Ивахненко и др. - Киев. : Наукова думка, 1975.

12. Кристофидес, Н. Теория графов. Алгоритмическая наука / Н. Кристофидес. - М. : Мир, 1978.

13. Кузин, Л. Т. Основы кибернетики / Л. Т. Кузин. - М. : Энергия, 1973.

14. Курицкий, Б. Я. Оптимизация вокруг нас / Б. Я. Курицкий. - Л. : Машиностроение, 1989.

15. Кушниренко, А. Г. Основы информатики и вычислительной техники / А. Г. Кушниренко. - М. : Просвещение, 1991.

16. Лященко, И. Н. Линейное и нелинейное программирование / И. Н. Ляшенко. - Киев. : Виша школа, 1975.

17. Мудров, А. Е. Численные методы для ПЭВМ на языках Бейсик,Фортран и Паскаль / А. Е. Мудров. - М. : МП “Раско”, 1992. 18. Никишев, В. К. Математический справочник-словарь по информатике и вычислительной технике / В. К. Никишев. - Чебоксары. : ЧГПИ, 1994.

19. Никишев, В. К. Использование IBM в учебном процессе / В. К. Никишев. - Чебоксары. : ЧГПИ, 1994.

20. Никишев, В. К. Основы программирования на ЦВМ и СЦВМ / В. К. Никишев. - Казань. : Министерство обороны, 1985.

21. Самарский, Н. Н. Компьютеры, модели, вычислительный эксперимент. Введение в информатику с помощью математического моделирования / Н. Н. Самарский. - М. : Наука, 1988.

22. Фомин, С. В. Математические проблемы в биологии. / С. В. Фомин. - М. : Наука, 1973.

23. Фридман, М. М. Наглядность и моделирование в обучении / М. М. Фридман. - М. : Знание, 1984.

24. Эберт, К. Компьютеры. Применение в химии / К. Эберт. - М. : Мир, 1988.


























Оглавление


1. Введение

2. Лабораторные работы по моделированию…………………4

2.1 Лабораторная работа №1

Тема: Методы исследования дифференциальных уравнений

с использованием разных информационных технологий…… 8

2.2 Лабораторная работа № 2

Тема: Свободное падение тел………………………………… 8

2.3 Лабораторная работа № 3

Тема: Исследование динамики объектов, брошенных

под углом к горизонту………………………………………….40

2.4 Лабораторная работа №4

Тема: Исследования динамики

популяций хищника и жертвы………………………………….68

2.5 Лабораторная работа №5

Тема: Моделирование движения

заряженных частиц и небесных тел……………………………81 2.6 Лабораторная работа № 6

Тема: Исследование динамики полета ракеты.……..………..98

2.7 Лабораторная работа № 6

Тема: Разработка информационной модели студента……..115

2.8 Лабораторная работа № 8

Тема: Исследование экономических систем ……..………….127

Тема: Исследование систем массового обслуживания … 143

2.9 Лабораторная работа № 9

Тема: Исследование оптимальных систем……………………150

2.10 Лабораторная работа №10

Тема. Моделирование объектов методом

пространства состояния………………………………… 159

Литература








Учебное издание

Лабораторный практикум

по моделированию



Учебно-методическое пособие Подписано в печать_27.02.07 Формат 60х80/16.

Бумага писчая. Печать оперативная. Усл. печ. л. 18.7

Тираж 50 экз. Заказ N .

ГОУ ВПО «Чувашский государственный

университет им. И. Н. Ульянова»

Отпечатано на участке оперативной полиграфии

ГОУ ВПО «Чувашский государственный

университет им. И.Н. Ульянова»

428000, Чебоксары,


















Автор
Дата добавления 08.04.2016
Раздел Информатика
Подраздел Конспекты
Просмотров1230
Номер материала ДБ-016936
Получить свидетельство о публикации

Выберите специальность, которую Вы хотите получить:

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

ПЕРЕЙТИ В КАТАЛОГ КУРСОВ

Похожие материалы

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