Математика \ Численные методы

Решение задачи Коши для системы обыкновенных дифференциальных уравнений, страница 2

y2[1] = opposed_method(y0[1], y0[0], 0.5*(Pi - (int)(Pi/dx) * dx), -0.05);

cout << "opposed_optimized: " << y2[0] << '\t' << y2[1] << endl;

//---------------------------------------------------------------------------------------------------------//Euler's method (for control):

y0[0] = 0.05*Pi;

y0[1] = 0.1*Pi;

dx = Pi/10000;

for(int i=0; i<=10000; i++)

{

y[0] = y0[0] + dx*0.1*y0[1];

y[1] = y0[1] - dx*0.05*y0[0];

y0[0] = y[0];

y0[1] = y[1];

}

y0[0] = 0.05*Pi;

y0[1] = 0.1*Pi;

dx=dx/2;

for(int i=0; i<=10000; i++)

{

y2[0] = y0[0] + dx*0.1*y0[1];

y2[1] = y0[1] - dx*0.05*y0[0];

y0[0] = y[0];

y0[1] = y[1];

}

if(abs(1.0*(y2[0]-y[0])/7) > e && abs(1.0*(y2[1]-y[1])/7) > e)

{

cout << "accuracy is not received..:( =->";

}

cout << "\nEuler's: " << y[0] << '\t' << y[1] << endl;

cout << "accuracy: y[0]:" << abs(1.0*(y2[0]-y[0])/7) << "  y[1]: " << abs(1.0*(y2[1]-y[1])/7) << endl;

//the thirdpart of the task : the step is variable :------------------------------------------------------------cout << "\nstep is variable:\n";

long double x = 0;

dx = 0.0001;//pow(e/delta, (long double)1.0/3);

y0[0] = 0.05*Pi;

y0[1] = 0.1*Pi;

long double r1,r2,r;

while(true)

{

//с шагом dx

y[0] = my_method(y0[0], y0[1], dx, 0.1);

y[1] = my_method(y0[1], y0[0], dx, -0.05);

//с шагом dx/2

y2[0] = my_method(y0[0], y0[1], 0.5*dx, 0.1);

y2[1] = my_method(y0[1], y0[0], 0.5*dx, -0.05);

y2[0] = my_method(y2[0], y2[1], 0.5*dx, 0.1);

y2[1] = my_method(y2[1], y2[0], 0.5*dx, -0.05);

//локальные погрешности для них:

r1 = 4.0*(y2[0]-y[0])/3;

r2 = 4.0*(y2[1]-y[1])/3;

r = max(r1, r2);

//оценка локальной погрешности:

if (r > 4*dlt)

{

dx = dx/2; // уменьшаем шаг в 2 раза и проходим этот шаг еще раз, т.к. точность не достаточна

}

else if ((r > dlt) && (r <= 4*dlt))

{              // уменьшаем шаг и делаем следующий

x += dx;

y0[0] = y2[0];

y0[1] = y2[1];

dx = dx/2;

}

else if ((r >= dlt/8) && (r < dlt))

{              // оставляем шаг и делаем следующий

x += dx;

y0[0] = y[0];

y0[1] = y[1];

}

else if (r < dlt/8)

{              // увеличиваем шаг в 2 раза и делаем следующий

x += dx;

y0[0] = y[0];

y0[1] = y[1];

dx = dx*2;

}

if(Pi-x < dx)

{

y[0] = my_method(y0[0], y0[1], Pi-x, 0.1);

y[1] = my_method(y0[1], y0[0], Pi-x, -0.05);

break;

}

if(x>=Pi) break;

}

cout << "1(mine): " << "y[0]=" << y[0] << '\t' << "y[1]=" << y[1] << endl;

//-----------opposed method:-------------------------x = 0;

dx = 0.0001;

y0[0] = 0.05*Pi;

y0[1] = 0.1*Pi;

r1=0; r2=0; r=0;

while(true)

{

//с шагом dx

y[0] = opposed_method(y0[0], y0[1], dx, 0.1);

y[1] = opposed_method(y0[1], y0[0], dx, -0.05);

//с шагом dx/2

y2[0] = opposed_method(y0[0], y0[1], 0.5*dx, 0.1);

y2[1] = opposed_method(y0[1], y0[0], 0.5*dx, -0.05);

y2[0] = opposed_method(y2[0], y2[1], 0.5*dx, 0.1);

y2[1] = opposed_method(y2[1], y2[0], 0.5*dx, -0.05);

//локальные погрешности для них:

r1 = 4.0*(y2[0]-y[0])/3;

r2 = 4.0*(y2[1]-y[1])/3;

r = max(r1, r2);

//оценка локальной погрешности:

if (r > 4*dlt)

{

dx = dx/2; // уменьшаем шаг в 2 раза и проходим этот шаг еще раз, т.к. точность не достаточна

}

else if ((r > dlt) && (r <= 4*dlt))

{              // уменьшаем шаг и делаем следующий

x += dx;

y0[0] = y2[0];

y0[1] = y2[1];

dx = dx/2;

}

else if ((r >= dlt/8) && (r < dlt))

{              // оставляем шаг и делаем следующий

x += dx;

y0[0] = y[0];

y0[1] = y[1];

}

else if (r < dlt/8)

{              // увеличиваем шаг в 2 раза и делаем следующий

x += dx;

y0[0] = y[0];

y0[1] = y[1];

dx = dx*2;

}

if(Pi-x < dx)

{

y[0] = opposed_method(y0[0], y0[1], Pi-x, 0.1);

y[1] = opposed_method(y0[1], y0[0], Pi-x, -0.05);

break;

}

if(x>Pi) break;

}

cout << "2(opposed): "<< "y[0]=" << y[0] << '\t' << "y[1]=" << y[1] << endl;

system("pause");

return 0;

}

файл defs.h:

#include <iostream>

#include <math.h>

using namespace std;

#ifndef Pi

#define Pi 3.141592653

#endif //Pi

#define e 0.0001

#define dlt 0.00001

long double my_method(long double y0, long double y1, long double step, long double par)

{

long double k1, k2;

k1 = step*par*y1;

k2 = step*par*(y1+0.1*k1);

return (y0 - 4.0*k1 + 5.0*k2);

}

long double opposed_method(long double y0, long double y1, long double step, long double par)

{

long double k1 = step*par*y1;

long double k2 = step*par*(y1 + 0.5*k1);

long double k3 = step*par*(y1 + 0.5*k2);

long double k4 = step*par*(y1 + k3);

return (y0 + 1.0/6*(k1 + 2*k2 + 2*k3 + k4));

}

long double max(long double x1, long double x2)

{

if (x1 > x2) return x1;

else return x2;

}

результатпрограммы (output):

-----my method----accuracy is not received..:( =->

mine: 0.2077863235      0.2992676418

accuracy: y[0]: 0.01614592185  y[1]: 0.006249335065

with optimized step:

mine_optimized: 0.2512873861    0.2819184967

-----opposed method----accuracy is not received..:(

opposed: 0.2054372733   0.2998917344

accuracy: y[0]: 0.003055854802  y[1]: 0.001197850612

with optimized step:

opposed_optimized: 0.2511189736 0.2819649253

Euler's: 0.2511155404   0.2819646966

accuracy: y[0]:6.327272995e-007  y[1]: 2.817509774e-007

step is variable:

1(mine): y[0]=0.2576983304      y[1]=0.2849461427

2(opposed): y[0]=0.2579070444   y[1]=0.2849267119

Для продолжения нажмите любую клавишу . . .

4. Использованная литература:

1. Олемской И.В. Методические указания по вычислительному практикуму (численные методы решения задачи Коши для СОДУ)

2.    Конспект лекций по курсу «Численные методы»

Скачать файл