Решение о.д.у. Задача Коши

Бесплатно!

Решение о.д.у. Задача Коши

Применяя метод Рунге-Кутта 3-го порядка точности и используя формулу 3.13  из пособия В.А. Юрко «Численные методы», найти решение следующей задачи  в точках 2.0, 4.0, 8.0. Нарисовать фазовый портрет на плоскости.

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

dr/dt=2r(1-r/500)-arf

df/dt=-f+arf   r(0)=r0, f(0)=f0, где t-время, r(t)-число кроликов,  f(t)-число лис и a- константа. Для a=0.01 найти решение при следующих               начальных условиях: r0=15, f0=22, t<=10.

Вычисления провести с точностью 1.

Алгоритм:

Формулы метода Рунге-Кутта:

Задаем начальные условия, определяем функцию правой части. Затем описываем функцию для вычисления начального шага, удовлетворяющего нашей заданной точности. Для этого задаем начальный шаг h, и для него считаем коэффициенты K1, K2 и находим w_1.  Далее берем шаг h1равный половине h и считаем для него ,получаем w_2. Сравниваем и смотрим удовлетворяет ли он заданной точности. Если да, то берем его, в противном случае делим шаг еще пополам. Найдя шаг и используя формулы метода Рунге-Кутта, вычисляем коэффициенты . Затем используя рекуррентную формулу для нахождения yj, находим значения y1 и y2 для соответствующих точек t.

#include
#include
#include
#include
#include
#include
using namespace std;

double a,h;
pair<double,double> f (double t, pair<double,double> y)
{
return make_pair(-y.first * y.second,y.first * y.first+y.second * y.second);
}

pair<double,double> euler (double t, pair<double,double> y)
{
pair<double,double> tmp(f(t,y));
return make_pair(tmp.first * h + y.first, tmp.second * h + y.second);
}

pair<double,double> f0 (double t, pair<double, double> y)
{
pair<double, double> a1(f(t, y));
return make_pair(h * a1.first, h * a1.second);
}
pair<double,double> f1 (double t, pair<double,double> y)
{
pair<double,double> a1(f0(t, y));
pair<double,double> a2(f(t+h/2,make_pair(y.first+a1.first/2,y.second+a1.second/2)));
return make_pair(a2.first * h, a2.second * h);
}
pair<double,double> f2 (double t, pair<double,double> y)
{
pair<double,double> a0(f0(t, y));
pair<double,double> a1(f1(t, y));
pair<double,double> a2(f(t+h,make_pair(y.first — a0.first + a1.first * 2, y.second — a0.second + a1.second * 2)));
return make_pair(a2.first * h, a2.second * h);
}
pair<double,double> rungekutta (double t, pair<double,double> y)
{
pair<double,double> a0(f0(t, y));
pair<double,double> a1(f1(t, y));
pair<double,double> a2(f2(t, y));
return make_pair(y.first + (a0.first + 4 * a1.first + a2.first)/6.0, y.second + (a0.second + 4 * a1.second + a2.second)/6.0);
}

double norm (const vector<pair<double,double> > &a, const vector<pair<double,double> > &b)
{
double res = 0;
for (int i = 0; i < a.size(); ++i) { double tmp = fabs(a[i].first — b[2 * i].first); if (tmp > res) res = tmp;
tmp = fabs(a[i].second — b[2 * i].second);
if (tmp > res) res=tmp;
}
return res;
}
vector<pair<double,double> > integr(int n, pair<double,double> nu)
{
h = 0.5/n;
vector<pair<double,double> > result;
result.push_back(nu);
for (int i = 0; i < n; ++i)
{
result.push_back(rungekutta(i * h,result[i]));
}
return result;
}

void printv(vector<pair<double, double> > a)
{
cout << endl;
for (int i = 0; i < a.size(); ++i)
cout << a[i].first << » » << a[i].second << endl;
cout << endl;
}

vector<pair<double,double> > runge(pair<double,double> nu, double eps)
{
int n = 10;
vector<pair<double, double> > v1(integr(n, nu));
while (true)
{
n *= 2;
vector<pair<double,double> > v2(v1);
v1 = integr(n, nu);
double tmp = norm(v2, v1);
if (tmp < eps)
{
return v1;
}
}
return vector<pair<double, double> >();
}
int main()
{
double eps = 0.01;
pair<double,double> nu(-2, 0);
vector<pair<double, double> > res(runge(nu, eps));
ofstream out(«output.txt»);
out << res[res.size()*0.1].first << «\t» << res[res.size()*0.1].second << endl;
out << res[res.size()*0.3].first << «\t» << res[res.size()*0.3].second << endl;
out << res[res.size()*0.4].first << «\t» << res[res.size()*0.4].second << endl;
}

Результат:

x1                  x2

-1.99002     0.200002

-1.91197     0.600488

-1.8461       0.802055

Детали:

Тип работы: Задачи, Контрольная работа

Предмет: Компьютерные науки

Год написания: 2011

Добавить комментарий

Ваш email не будет показан.

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