исследование моделей взаимодействия двух видов "Хищник-жертва" и "Конкуренция" с использованием метода Рунге -Кугга 4 порядка и решение задачи Коши.

Исследование полученных результатов подтвердило более высокую точность метода Рунге-Кутта 4 порядка по сравнению с методом Эйлера.

Рис.6.1.Внешний вид интерфейса программы
Список литературы
1. Базыкин А.Д. Математическая биофизика взаимодействующих популяций. - М.: Наука, 1985. – 184 с.
2. Ортега Дж., Пул У. Введение в численные методы решения дифференциальных уравнений/ Пер. с англ.; Под ред. А.А. Абрамова – М.: Наука. Гл. ред. физ.-мат. лит., 1986. – 288 с., ил.
3. Вольтерра В. Математическая теория борьбы за существование. – М.: Наука, 1976. – 288 с.
4. Алгоритмы и программы проектирования автоматических систем/ П.Д. Крутько, А.И. Максимов, Л.М. Скворцов; Под ред. П.Д. Крутько. – М.: Радио и связь, 1988. – 306 с., ил.
Приложение
void
nextRK4(eco *s,double h) // Реализация метода Рунге-Кутта 4-го порядка
{
int d=s->x.getd();
double t=s->time;
column k1(d),k2(d),k3(d),k4(d);
column x1;

x1=s->x;
k1=s->fr(t,x1)*h;
t += h/2.;
k2=s->fr(t,x1+k1/2.)*h;
k3=s->fr(t,x1+k2/2.)*h;
t += h/2;
s->x=x1+(k1+k2*2.+k3*2+s->fr(t,x1+k3)*h)/6.;
s->time += h;
}

void
nextEil(eco *s,double h) // Реализация метода Эйлера
{
s->x=s->x+s->fr(s->time,s->x)*h;
s->time += h;
}


void CCgrafDlg::OnButton1()
{
// TODO: Add your control notification handler code here
// процедура интегрирования
double hp; // шаг фиксации результата
long icount; // количество точек интегрирования
int pcount; // количество точек для вывода
double *t; // время в зафиксированных точках
double *result; // переменные состояния в зафиксированных точках
double min,max; // пределы изменения переменных
double eps,e1,e2; // оценка погрешности
char buf[20];
int i,j;

icount=(int)(itime/h)+1; // определение шага фиксации результатов
if(icount < 10001){ // берем не более 10000 точек
pcount=icount;
hp=h;
}
else{
pcount=10001;
hp=itime/10000.;
}

t= new double[pcount];
result= new double[pcount*2];

m_Pr.SetPos(0); // инициализация прогрессора
m_Pr.SetRange(0,pcount);

column x0(2);
x0.set(0,x[0]); // начальные условия
x0.set(1,x[1]);

eco s(x0);
eco s2(x0);
// цикл интегрирования
eps=0.;
for(i=0,j=0; i < icount; i++){
double mi=s.getx(0) < s.getx(1) ? s.getx(0) : s.getx(1);// попутно определяем минимум
double ma=s.getx(0) > s.getx(1) ? s.getx(0) : s.getx(1);// и максимум для построения
if(i == 0){ // первую точку включаем безусловно
min=mi;
max=ma;
t[j]=s.gettime();
result[j]=s.getx(0);
result[pcount+j]=s.getx(1);
j++;
}
else{
min=min < mi ? min : mi; // попутно определяем минимум
max=max > ma ? max : ma; // и максимум для построения
}
if(s.gettime()+h >= t[j-1]+hp && j < pcount-1){// включаем в результаты очередную
t[j]=s.gettime()+h; // точку
result[j]=s.getx(0);
result[pcount+j]=s.getx(1);
j++;
m_Pr.SetPos(j); // двигаем прогрессор
}
if(iMeth){ // выбран метод Эйлера
nextEil(&s,h);
nextEil(&s2,h/2); // интегрируем кроме того с шагом h/2
nextEil(&s2,h/2);
e1=fabs(s.getx(0)-s2.getx(0)); // чтобы определить погрешность
e2=fabs(s.getx(1)-s2.getx(1));
eps=eps > e1 ? eps > e2 ? eps : e2 : e1;
}
else{ // выбран метод Рунге-Кутта
nextRK4(&s,h);
nextRK4(&s2,h/2); // интегрируем кроме того с шагом h/2
nextRK4(&s2,h/2);
e1=fabs(s.getx(0)-s2.getx(0)); // чтобы определить погрешность
e2=fabs(s.getx(1)-s2.getx(1));
eps=eps > e1 ? eps > e2 ? eps : e2 : e1;
}
}
if(!iMeth)
eps /= 15.; // для метода 4-го порядка
sprintf(buf,"%11.5g",eps); // выводим полученную погрешность
m_Eps.SetWindowText(buf);
for(;j < pcount; j++){ // если остались незаполненные результаты, t[j]=s.gettime(); // тиражируем последние значения

result[j]=s.getx(0);
result[pcount+j]=s.getx(1);
m_Pr.SetPos(j); // доводим прогрессор
}
if(min > 0.) // догоняем минимумы и максимумы до целых
min=0.;
else
min=floor(1.1*min);
max=ceil(1.1*max);
m_Chart2d.AllocSerie(pcount); // выводим график
m_Chart2d.SetRange(0,itime,min,max);
for(i=0; i < pcount; i++){
m_Chart2d.SetXYValue(t[i],result[i],i,0);
m_Chart2d.SetXYValue(t[i],result[pcount+i],i,1);
}
m_Chart2dh.AllocSerie(pcount); // выводим фазовый портрет
m_Chart2dh.SetRange(min,max,min,max);
for(i=0; i < pcount; i++)
m_Chart2dh.SetXYValue(result[i],result[pcount+i],i,0);
m_Chart2dh.Invalidate();
m_Chart2d.Invalidate();
UpdateData(FALSE);

delete[] t; // освобождаем память
delete[] result;
}















18



y

A

e

B

y=y(x)

yi+1

h

yi

x

O

хi

xi+1

α