Смотреть ещё
682
методические разработки по информатике
Перейти в каталогМИНОБРНАУКИ РОССИИ
Федеральное государственное бюджетное образовательное
учреждение высшего профессионального образования
«Чувашский государственный университет имени И. Н. Ульянова»
В. К. Никишев
Математическое
моделирование
Лабораторный практикум
Чебоксары
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/
Блок-схема алгоритма
Отчет по лабораторной работе
Рис.1. Форма 2. Автор Рис. 2. Форма 4. Условие
Рис. 3. Форма 1. Титульный лист
Форма для исследования объекта
Рис. 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;
Рис 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);
}
while (x<=b);
}};}};}
Рис 6. Форма для автора.
Рис 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
Рис 8. Титульный лист Рис. 9. Условия задачи
Рис 10. Автор программы
Рис. 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++
Рис. 12. Титульный лист
Рис 13. Автор программы Рис 14. Условие задачи
Лист 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.
Рис. 15. Титульный лист и условие задачи
Рис. 16. Задачи моделирования
Рис. 17. Результаты моделирования
Блок-схема алгоритма
Программа на 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)
Рис. 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);
Рис. 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);
Рис. 21. Результаты исследования
Пример 4
Условия задачи. Исследовать объект, динамика которого описывается дифференциальным уравнением
ӱ + 0.2*ý+ 4*y=2
методом визуального моделирования
Исходное уравнение преобразуем к виду
Ӱ=2-0.2*ý-4*y
Рис. 22. Визуальная схема моделирования
Рис. 23. Результаты моделирования
Рис. 24. Визуальная схема моделирования
Рис. 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
Рис. 25. Результаты моделирования
Пример
Исследование в среде MatLab ( Simulink).
Условие задачи: составить алгоритм и проект:
моделирования объекта, динамика которого описывается
дифференциальным уравнением 2 порядка
y’’+axy’-y=0,4.
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.
Модель объекта
Параметры k1 и k2 зависят от свойств среды и геометрии тела. Например, для шарика k1 =6μπr (формула Стокса), где μ– динамическая вязкость среды, r– радиус шарика. Так, для воздуха μ = 0,0182 Н*с*м-2 при t=20o C и давлении 1 атм. Соответственно для воды μ = 1,002 Н*с*м-2. Для малых скоростей составляющей k1 можно пренебречь и Fс = k2*v2, где k2 = 0,5сsρ ; ρ – плотность среды ; с– коэффициент лобового сопротивления, который зависит от формы тела. Так, для диска в форме прямоугольника с=1,11, для сферы с=1,33 и с=0,55 в зависимости от направления расположения сферы, для шара с=0,4, а для каплевидных тел с=0,045 с=0,1 в зависимости от направления тела.
.
Пример. Исследовать падение шарика радиуса 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 - коэффициенты внешней среды.
2. - площадь сечения тела, поперечного по отношению к потоку
3. плотность шарика;
4. – плотность воздуха.
5.
6. c = 0.4 – коэффициент лобового сопротивления шарика;
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;
form5.hide;
form1.Show;
end;
Пример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
Исследование системы с использованием программы MatLab.
Графики зависимостей v(t), h(t) для r=5mm c
помощью решателей ode15,ode45
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<R<5mm
Simulink
Пример3. Исследовать движение исследовательского зонда вертикально вверх с летящего самолета
Условие
Промоделировать движение исследовательского зонда, «выстреленного» вертикально вверх с летящего над землей самолета. В верхней точке траектории над зондом раскрывается парашют, и он плавно спускается на землю.
Входные параметры.
Математическая модель свободного падения тела - уравнение второго закона Ньютона с учетом двух сил, действующих на тело -силы тяжести и силы сопротивления среды. Движение является одномерным; проецируя силу, скорость и перемещение на ось, направленную вертикально вверх, получаем
Входные параметры модели:
¦ начальная высота тела;
¦ начальная скорость тела;
¦ величины, определяющие коэффициенты сопротивления среды
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]);
Пример. Исследование динамики объектов, брошенных под углом к горизонту.
Цель занятия:
Получить практические навыки исследования объектов, динамика которых описывается дифференциальными уравнениями 2-го порядка с использованием языков программирования :Delphi, VC++, VC# и информационных технологий MatLab.
Задачи занятия:
1. Разработка алгоритмов
2. Построение семейства кривых y(x), dy/dx(x) при параметрах a-const, a-var
3. Получение зависимостей y(a), dy/da при x-const
4. Анализ результатов исследований.
Задания для проведения лабораторной работы
Динамика объекта описывается системой дифференциальных уравнений
mdvx/dt = -Fccosq; =-Fcvx/v0
mdvy/dt = -m*g - Fcsinq,=-m*g - Fcvy/v0,
где Fc =k1V+ k2*V2 - линейная и квадратичная составляющая сопротивления воздуха.
2. Динамика объекта может представляется и алгебраическими уравнениями вида
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<y0) then begin
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<y0) then begin
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<shape17.Top+shape17.height)or
(y>shape16.Top+shape16.height)) then break;
end;
if x>=shape22.Left then begin
if (y<shape22.Top)or
(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.left:=paintbox2.left+x;
paintbox3.left:=paintbox3.left+x;
paintbox4.left:=paintbox4.left+x;
shape15.Left:=shape15.Left+x;
shape16.Left:=shape16.Left+x;
shape17.Left:=shape17.Left+x;
shape18.Left:=shape18.Left+x;
shape19.Left:=shape19.Left+x;
shape20.Left:=shape20.Left+x;
shape23.Left:=shape23.Left+x;
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.left:=paintbox3.left+x;
shape17.width:=shape15.width+20;
shape17.left:=shape15.left-10;
paintbox4.width:=shape15.width+20;
paintbox4.left:=shape15.left-10;
shape18.width:=shape17.width div 5;
shape19.width:=shape17.width div 5;
shape18.left:=shape17.left;
shape19.left:=shape17.left+shape17.width div 5*2;
shape20.left:=shape17.left+shape17.width div 5*4;
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.top:=shape17.top+y;
shape18.top:=shape18.top+y;
shape19.top:=shape19.top+y;
shape16.top:=shape16.top+y;
shape20.top:=shape20.top+y;
shape15.top:=shape15.top+y;
paintbox4.top:=paintbox4.top+y;
paintbox3.top:=paintbox3.top+y;
paintbox2.top:=paintbox2.top+y;
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.top:=shape22.top+y;
paintbox5.top:=paintbox5.top+y;
paintbox6.top:=paintbox6.top+y;
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.top:=paintbox6.top+y;
end;
end;
procedure TForm1.PaintBox1MouseUp(Sender: TObject; Button: TMouseButton;
Shift: TShiftState; X, Y: Integer);
begin
dr:=false;
end;
end.
Пример модели движения небесных тел
По закону всемирного тяготения сила притяжения, действующая между двумя телами, пропорциональна их массам и обратно пропорциональна квадрату расстояния между ними.
В этой формуле G=6,67*10-11 m3/(kг*с2) -гравитационная постоянная. В данной системе координат начало координат расположено на одном теле M. Уравнения, описывающие движение тела m в указанной системе координат, имеет вид
В проекциях на оси координат уравнение можно привести с следующей системе из 4-х дифференциальных уравнений
Пример 2. Исследовать выполнение второго закона Кеплера при движении различных планет.
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.
2.4 Лабораторные работы по разработке имитационных моделей
Пример. Разработка информационной модели студента ( учащихся)
Цель занятия: Получить практические навыки в разработке информационных систем в среду БД Access , Excel
Задачи занятия:
2. Создать графическую модель студента с отображением диаграммы оценок по темам конкретной дисциплины
Информационная модель студента
Предметная информационная модель знаний студентов курса или факультета позволит получить следующие характеристики учебного процесса:
- качество обучения в целом по университету, факультету, курсу и успеваемость отдельных студентов;
-стабильность обучения по отдельным предметам, темам и дисциплинам;
- прогнозирование результатов обучения в предстоящем семестре;
-систематическое уточнение модели знаний студентов с использованием тестовых программ;
- разработка индивидуальных вопросов и задач на основе моделей знаний студентов;
- сравнительные оценки состояния учебного процесса как на факультете, так и по университету в целом.
Пример имитационной модели знаний студентов в среде Excel
Пример имитационной модели учащихся школы в среде Access
2.5 Разработка моделей транспортных задач
Пример. «Размещения предприятий»
Исходные условия. Пусть в нескольких пунктах расположены предприятия, производящие некоторый продукт. В пунктах размещаются потребители готовой продукции с соответствующими потребностями . Затраты на производство единицы продукта в пункте равны , обьем производства в этом пункте равен , а затраты по транспортировке единицы продукта из i в j равны . Количество перевезенных продуктов из i в j составляет единиц. Уравнения модели:
1.(перевозится неотрицательное количество продукта);
2. (выпускаемое количество продукта не
больше возможного обьема производства и равно вывозимому количеству продукта);
3. , j=1,...,n (потребности потребителей),
(каждый пункт потребления получает столько, сколько ему требуется).
Конечная цель. Необходимо получить минимальные затраты на производство и на транспортировку продукции, т.е. обеспечить минимум функции
,
где - затраты на производство,
- затраты на транспортировку.
Целевая функция состоит из аналогичной целевой задачи. Однако, различие состоит в том, что здесь добавляются затраты на производство. Для однотипности подхода к решению этих двух моделей добавляется фиктивный пункт с характеристикой: потребность равна разности между возможным обьемом производства продукта и суммарной потребнос
,
где - затраты на производство,
- затраты на транспортировку.
Целевая функция состоит из аналогичной целевой задачи. Однако, различие состоит в том, что здесь добавляются затраты на производство. Для однотипности подхода к решению этих двух моделей добавляется фиктивный пункт с характеристикой: потребность равна разности между возможным обьемом производства продукта и суммарной с
Пример Моделирование системы планирования
на основе метода сетевого графа
Задачи исследования
Необходимо спланировать временный график работ в соответствии с графом и вычислить критический путь выполнения всех работ
Исходные данные
Задание: Дана схема, нужно найти критический путь
Пример. Планирование производства товаров на основе модели получения максимальной прибыли с использованием метода линейного программирования
Условие. Пусть предприятие производит столы и стулья. Расход ресурсов на их производство и прибыль от их реализации представлены в таблице. столы стулья объем ресурсов
Расход древесины
на изделие, м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 ( целевая функция)
unit
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. Анализ результатов исследований с вычислением параметров и эффективности систем.
Задания для исследований в разделе Задания
Пример 1. Исследовать параметры обслуживания покупателей в магазине с одним продавцом
Программа на Делфи
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.
2.9
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<y Then begin {если еще есть пришедшие}
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.Left:=Image1.Left-30;
If Image1.Left<112 Then begin
Timer2.Enabled:=False;
Image1.Visible:=False;
Image1.Left:=208;
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<kon do
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. Анализ результатов исследований с вычислением параметров
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<rast then begin
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 then begin
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.
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.
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
Математическая модель объекта.
Исходное дифференциальное уравнение имеет вид
Представим исходное уравнение в нормальной форме Коши, вводя обозначения
В результате получим систему уравнений первого порядка
С целью определения матриц уравнения параметров состояния преобразуем систему:
матрицы коэффициентов соответственно равны
Уравнение параметров состояния будет иметь вид
Решение данных примеров в системе 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
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
Пример. Для дифференциального уравнения
С целью определения матриц уравнения параметров состояния преобразуем систему:
Пример 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 //сетка
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<L<1.
Более строго соотношение может быть получено из системы дифференциальных уравнений:
Модель “хищник - жертва”
Имеются популяции двух видов, которые представляются отношениями
xi+1 = xi + a1xi - b1xi2 - g1xiyi , x0 = c1,
yi+1 = yi - a2yi + b2yi2 - g2xiyi , y0 = c2,
где xi - численность (плотность) жертв,
yi - численность хищников,
g1 - коэффициент защиты жертв,
g2 - коэффициент прожорливости хищников.
Рост опухоли
Раковая опухоль обычно увеличивается экспоненциально в соответствии с дифференциальным уравнением:
, где v - размер опухоли,
с,a,b - константы.
Определить, при каких значениях параметров С существует предельный размер опухоли. Выяснить, при каких значениях С рост опухоли не превосходит некоторой конечной величины.
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. Значения остальных параметров фиксировать по усмотрению.
23.Промоделировать процесс распространения инфекции стригущего лишая по участку кожи размером 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. Кратчайший маршрут
Две магистрали пересекаются под острым углом, внутри которого расположены два населенных пункта. Как проложить кратчайший маршрут автобуса, соединяющий два данных пункта и имеющий выезды к каждой из двух магистралей в заданном порядке?
Задания для метода « Линейное программирование»
Индивидуальные задания
- модель движения материальной точки Аристотеля и Ньютона;
- модель Солнечной системы Птолемея, Коперника, Кеплера;
- простейшую демографическую модель;
- модель многоотраслевой экономики Леонтьева;
- модель процесса распространения эпидемий;
- модель динамики численности биологических популяций;
- модель относительных движений в классической механике;
- модель остывания нагретых тел в атмосфере;
- математическую модель процесса загрязнения воды;
- модель колебательных процессов в физике.
- модель биологической системы «хищник-жертва»;
- модель биологической системы конкурирующих популяций;
- модель поведение динамики многочастичной системы;
- модели марковских случайных процессов;
- модель поведение динамической системы, описываемой
уравнениями Колмогорова.;
-модель управления различными летательными объектами;
- модель исследования систем на основе матричных методов пространства состояния;
- модели экспертных систем;
- модели управления педагогическими системами;
- модели геоинформационных систем .
ЛИТЕРАТУРА
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, Чебоксары,
В нашем каталоге доступно 74 699 рабочих листов
Перейти в каталогПолучите новую специальность за 2 месяца
Получите профессию
за 6 месяцев
Пройти курс
Рабочие листы
к вашим урокам
Скачать
6 665 176 материалов в базе
Настоящий материал опубликован пользователем Никишов Вячеслав Константинович. Инфоурок является информационным посредником и предоставляет пользователям возможность размещать на сайте методические материалы. Всю ответственность за опубликованные материалы, содержащиеся в них сведения, а также за соблюдение авторских прав несут пользователи, загрузившие материал на сайт
Если Вы считаете, что материал нарушает авторские права либо по каким-то другим причинам должен быть удален с сайта, Вы можете оставить жалобу на материал.
Удалить материалВаша скидка на курсы
40%Курс профессиональной переподготовки
500/1000 ч.
Курс повышения квалификации
72 ч. — 180 ч.
Курс повышения квалификации
72 ч. — 180 ч.
Курс профессиональной переподготовки
300 ч. — 1200 ч.
Мини-курс
4 ч.
Мини-курс
5 ч.
Оставьте свой комментарий
Авторизуйтесь, чтобы задавать вопросы.