Код:
#include <string>
#include <iostream>
#include <cstring>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <time.h>
#include <conio.h>
#include <stdlib.h>
using namespace std;
struct Var // переменные системы первого порядка
{
double x, v;
};
/*// арифметические операции над структурами переменных:
Var operator+ (Var& A, Var& B);
Var operator* (double k, Var& A);*/
struct Par // структура параметров задачи
{
double omega0, gamma, F, omega;
};
/*// функция, задающая внешнюю силу
double Ext_force(Par& P, double t);
// функция, вычисляющая правые части системы первого порядка:
Var RHS(Var& V, Par& P, double t);
// шаг метода трапеций:
void Trapezium_step(Var& V, Par& P, double& t, double st);
// функции, вычисляющие точное решение для используемого вида внешней силы:
double Exact_x(Var& V0, Par& P, double t);
double Homogeneous_x(Var& V0, Par& P, double t);
double Partial_0_x(Par& P, double t);
*/
Var operator+ (Var& A, Var& B)
{ // сложение двух структур переменных
Var Res;
Res.x = A.x + B.x;
Res.v = A.v + B.v;
return Res;
}
Var operator* (double k, Var& A)
{ // умножение структуры переменных на число
Var Res;
Res.x = k * A.x;
Res.v = k * A.v;
return Res;
}
Var RHS(Var& V, Par& P, double t);
double Homogeneous_x(Var& V0, Par& P, double t);
double Partial_0_x(Par& P, double t);
void Trapezium_step (Var &V, Par &P, double &t, double st)
{
Var V0 = V; // структура переменных на предыдущем шаге
Var F0 = RHS(V0, P, t); // значение правых частей на предыдущем шаге
Var V_prev = V;
// вспомогательная структура, хранящая значение вычисляемой новой структуры переменных, полученное на предыдущей итерации
do // цикл решения алгебраического уравнения
{
// методом простых итераций
V_prev = V;
V = V0 + 0.5 * st * ( F0 + RHS(V_prev, P, t + st) );
}
while ( fabsl(V.x - V_prev.x)+ fabsl(V.v - V_prev.v) > 0.00000001 );
// сравнение результатов текущей и предыдущей итераций в норме, равной сумме модулей компонент; вычисление останавливается, если норма разности меньше заданного значения
t += st; // сдвиг текущего времени на шаг
}
double Ext_force(Par& P, double t)
{
// внешняя сила
return P.F * cosl(P.omega * t);
}
Var RHS(Var& V, Par& P, double t)
{
// правые части системы первого порядка
Var Res;
Res.x = V.v;
Res.v = - P.gamma * V.v - pow(P.omega0, 2) * V.x + Ext_force(P, t);
return Res;
}
double Exact_x(Var& V0, Par& P, double t)
{
// частное решение неоднородного уравнения с заданными начальными условиями
return Homogeneous_x(V0, P, t) + Partial_0_x(P, t);
}
double Homogeneous_x(Var& V0, Par& P, double t)
{
// частное решение однородного уравнения с заданными начальными условиями
double Omega = sqrtl( pow(P.omega0, 2) - pow(P.gamma, 2) / 4. );
return expl(- P.gamma * t / 2.) * ( V0.x * cosl(Omega * t) + ( (V0.v + P.gamma * V0.x / 2.) / Omega ) * sinl(Omega * t) );
}
double Partial_0_x(Par& P, double t)
{
// частное решение неоднородного уравнения с нулевыми начальными условиями
double Omega = sqrtl( pow(P.omega0, 2) - pow(P.gamma, 2) / 4. );
double diff = pow(P.omega0,2) - pow(P.omega,2);
double exponent = expl(- P.gamma * t / 2.);
double factor = ( pow(P.omega0,2) + pow(P.omega,2) ) / (2. * P.omega * Omega);
return ( P.F / ( pow(diff, 2) + pow(P.gamma, 2) * pow(P.omega, 2) ) ) *( diff * ( cosl(P.omega * t) - exponent * cosl(Omega * t) ) + P.gamma * P.omega * ( sinl(P.omega * t) - factor * exponent * sinl(Omega * t) ) );
}
int main()
{
Var V;
V.x = 0.;
V.v = 0.;
Var V0 = V;
Par P;
P.omega0 = 1.;
P.gamma = 0.05;
P.F = 1.;
P.omega = 1.2;
double t = 0.;
double st = 0.001;
double x_ex = 0.; // точное решение
double max_diff = 0.;
// максимальное отклонение численного решения от точного (на данный момент времени)
ofstream x_res_file("x.res");
// файл численной траектории
ofstream x_exact_res_file("x_exact.res");
// файл аналитической траектории
double t1 = st;
double error_out_interval = 1.;
// вспомогательные переменные для вывода на экран значения ошибки с данным интервалом времени error_out_interval
int i = 0;
int n_steps_out = 100;
// вспомогательные переменные для вывода в файл значений x через каждые n_steps_out шагов
while (t < 200) // рассчитываем траекторию
{
// на интервале времени (0, 200)
Trapezium_step(V, P, t, st);
// шаг интегрирования
x_ex = Exact_x(V0, P, t);
// вычислить точное решение
i++;
if ( i == n_steps_out )
{
// вывод численного и точного решений в файлы
x_res_file << V.x << " ";
x_exact_res_file << x_ex << " ";
i = 0;
}
// подсчет отклонения от точного решения
if ( fabsl(V.x - x_ex) > max_diff)
max_diff = fabsl(V.x - x_ex);
if (t1 >= error_out_interval)
{
// вывод отклонения
cout << "\ntime: " << t
<< " error: " << max_diff << " ";
t1 = st;
}
else
t1 += st;
}
x_res_file.close();
x_exact_res_file.close();
return 0;
}
27.03.2016, 22:40
ошибка при запуске ;(
Похожие темы
