Исследование моделей взаимодействия двух видов «Хищник-жертва» и «Конкуренция» с использованием метода Рунге — Кугга 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+k½.)*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 < 10 001){ // берем не более 10 000 точек

pcount=icount;

hp=h;

}

else{

pcount=10 001;

hp=itime/10 000.;

}

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;

}

y

A

e

B

y=y (x)

yi+1

h

yi

x

O

хi

xi+1

α