我想要达到的结果
// -----------------------------------------------------------------
// main.cpp
// -----------------------------------------------------------------
// This is the main file for a program which trace the Basins of
// attraction for the piezo-magneto-elastic beam
// -----------------------------------------------------------------
// programmer:
// João Peterson (ligier.peterson@gmail.com)
// Americo Cunha (americo.cunha@uerj.br)
//
// last update: Jan, 2021
// -----------------------------------------------------------------
#include
#include
#include
#include // ofstream
#include // round
#include "header.h"
using namespace std;
void task_t2(double *d_params, int *i_params, double *param, double** AT, double** IC);
void task_t3(double *d_params, int *i_params, double *param, double** AT, double** IC);
void task_t4(double *d_params, int *i_params, double *param, double** AT, double** IC);
int main()
{
cout << "// -------------------- //\n";
cout << "// Basins of attraction\n";
cout << "// -------------------- //\n" << endl;
cout << "Please, wait...\n" << endl;
// ------------------------------------------- //
// Simulation parameters
// ------------------------------------------- //
double t0 = 0.0; // initial time
double t1 = 5e3; // final time
double h = 5e-2; // time step
long int Ndt = round((t1-t0)/h); // number of time steps
long int Nss = round(0.9*Ndt); // Steady-state
double tol = 0.01; // tolerance for comparison test
int Nat = 10; // number of attractors
double x0_min = -3.0; // min. of initial displacement
double x0_max = 3.0; // max. of initial displacement
double xdot0_min = -3.0; // min. of initial velocity
double xdot0_max = 3.0; // max. of initial velocity
int N_x = 50; // number of points
int N_xdot = 50; // number of points
double *x0 = new double[N_x]; // vector of x0 values
linspace(x0_min, x0_max, N_x, x0);
double *xdot0 = new double[N_xdot]; // vector of xdot0 values
linspace(xdot0_min, xdot0_max, N_xdot, xdot0);
double v0 = 0.0;
double ic[3];
ic[2] = v0;
double d_params [] = {t0, t1, h, tol, x0_min, x0_max, xdot0_min, xdot0_max, v0};
int i_params[] = {Ndt, Nss, Nat, N_x, N_xdot};
// ------------------------------------------- //
// Physical parameters
// ------------------------------------------- //
double ksi, chi, f, omega, lambda, kappa;
//double ksi, chi, f, lambda, kappa;
ksi = 0.01; // mechanical damping ratio
chi = 0.05; // dimensionless piezoelectric coupling term (mechanical)
f = 0.083; // dimensionless excitation amplitude
omega = 0.5; // dimensionless excitation frequency
lambda = 0.05; // dimensionless time constant reciprocal
kappa = 0.5; // dimensionless piezoelectric coupling term (electrical)
double param[] ={ksi, chi, f, omega, lambda, kappa};
// ------------------------------------------- //
// Allocating memory
// ------------------------------------------- //
//attractors output matrix
/*
This will create a matrix with dimensions [Ndt-Nss+1][3*Nat],
the data is meant to be stored in that way:
-- first attractor -- -- second attractor --
x0[0] xdot0[0] v0[0] x0[0] xdot0[0] v0[0] ...
x0[1] xdot0[1] v0[1] x0[1] xdot0[1] v0[1] ...
... ... ... ... ... ...
*/
double **AT = new double*[Ndt-Nss+1];
for(int i=0; i[i] = new double[3*Nat];}
for(int i=0; i[i][j] = 0;}
}
// IC output matrix
/*
This matrix is organized as follows:
dimensions: [N_x][N_xdot]
z_[0][0] z_[0][1] z_[0][2]
z_[1][0] z_[1][1] z_[1][2]
... ... ...
*/
double **IC = new double *[N_x];
for(int i=0; i[i] = new double[N_xdot];}
// response matrix (dimensions: [Ndt][3])
double **y = new double*[Ndt];
for(int i=0; i[i] = new double[3];}
double *Qvolt = new double[(Ndt-Nss)/50];
// ------------------------------------------- //
// Integrating...
// ------------------------------------------- //
thread t2(task_t2, d_params, i_params, param, ref(AT), ref(IC));
thread t3(task_t3, d_params, i_params, param, ref(AT), ref(IC));
thread t4(task_t4, d_params, i_params, param, ref(AT), ref(IC));
for(int nx=0; nx[0] = x0[nx]; // updating displacement
for(int nxdot=0; nxdot[1] = xdot0[nxdot]; // updating velocity
RKF4(ic, 3, param, t0, t1, h, pvi, y); // Integrate the dynamical system with RKF4
double **yss = new double *[Ndt-Nss+1]; // steady-state response
for(int i=0; i[i] = new double[3];
yss[i][0] = y[Nss+i-1][0];
yss[i][1] = y[Nss+i-1][1];
yss[i][2] = y[Nss+i-1][2];
}
for(int i=0; i<(Ndt-Nss)/50; i++){
Qvolt[i] = yss[(i*50)][2];
}
if(abs(z1test(Qvolt, (Ndt-Nss)/50)) < 0.2) // if regular
{
for(int nat=0; nat[0][nat*3] == 0) // if attractor is empty
{
fill_AT(yss, 0, Ndt-Nss, AT, nat*3, nat*3+2); // fill attractor
IC[nx][nxdot] = nat+1; // fill IC matrix
break;
}
else if(test(yss, Ndt-Nss+1, AT, nat*3, tol)) // if same attractor
{
IC[nx][nxdot] = nat+1; // fill IC matrix
break;
}
}
}
else // chaotic
{
IC[nx][nxdot] = 0; // fill IC matrix
}
for(int i=0; i[] yss[i];}
delete [] yss;
yss = 0;
}
}
// Join threads
t2.join();
t3.join();
t4.join();
// save information
ofstream ATfile("AT_[" + to_string(int(x0_min)) + "," + to_string(int(x0_max)) + "]_" +
"[" + to_string(int(xdot0_min)) + "," + to_string(int(xdot0_max)) + "]_" +
to_string(N_x) + "X" + to_string(N_xdot) + "_f" + to_string(int(f*1000)) +
"_O" + to_string(int(omega*100)) + ".dat", ios::trunc);
for(int i=0; i[i][j] << " ";
}
ATfile << endl;
}
ATfile.close();
cout << "ok";
ofstream ICfile ("IC_[" + to_string(int(x0_min)) + "," + to_string(int(x0_max)) + "]_" +
"[" + to_string(int(xdot0_min)) + "," + to_string(int(xdot0_max)) + "]_" +
to_string(N_x) + "X" + to_string(N_xdot) + "_f" + to_string(int(f*1000)) +
"_O" + to_string(int(omega*100)) + ".dat", ios::trunc);
for(int j=0; j[i][N_xdot-1-j] << " ";
}
ICfile << endl;
}
ICfile.close();
delete [] x0;
x0 = 0;
delete [] xdot0;
xdot0 = 0;
for(int i=0; i[] y[i];}
delete [] y;
y = 0;
for(int i=0; i[] IC[i];}
delete [] IC;
IC = 0;
for(int i=0; i[] AT[i];}
delete [] AT;
AT = 0;
delete [] Qvolt;
Qvolt = 0;
return 0;
}
void task_t2(double *d_params, int *i_params, double *param, double** AT, double** IC){
double t0 = d_params[0];
double t1 = d_params[1];
double h = d_params[2];
double tol = d_params[3];
double x0_min = d_params[4];
double x0_max = d_params[5];
double xdot0_min = d_params[6];
double xdot0_max = d_params[7];
double v0 = d_params[8];
int Ndt = i_params[0];
int Nss = i_params[1];
int Nat = i_params[2];
int N_x = i_params[3];
int N_xdot = i_params[4];
// vectors
double *x0 = new double[N_x]; // vector of x0 values
linspace(x0_min, x0_max, N_x, x0);
double *xdot0 = new double[N_xdot]; // vector of xdot0 values
linspace(xdot0_min, xdot0_max, N_xdot, xdot0);
double ic[3];
ic[2] = v0;
// response matrix (dimensions: [Ndt][3])
double **y = new double*[Ndt];
for(int i=0; i[i] = new double[3];}
double *Qvolt = new double[(Ndt-Nss)/50];
// calculating
for(int nx=1; nx[0] = x0[nx]; // updating displacement
for(int nxdot=0; nxdot[1] = xdot0[nxdot]; // updating velocity
RKF4(ic, 3, param, t0, t1, h, pvi, y); // Integrate the dynamical system with RKF4
double **yss = new double *[Ndt-Nss+1]; // steady-state response
for(int i=0; i[i] = new double[3];
yss[i][0] = y[Nss+i-1][0];
yss[i][1] = y[Nss+i-1][1];
yss[i][2] = y[Nss+i-1][2];
}
for(int i=0; i<(Ndt-Nss)/50; i++){
Qvolt[i] = yss[(i*50)][2];
}
if(abs(z1test(Qvolt, (Ndt-Nss)/50)) < 0.2) // if regular
{
for(int nat=0; nat[0][nat*3] == 0) // if attractor is empty
{
fill_AT(yss, 0, Ndt-Nss, AT, nat*3, nat*3+2); // fill attractor
//IC[nx*N_xdot+nxdot][0] = ic[0]; // fill IC matrix
//IC[nx*N_xdot+nxdot][1] = ic[1];
//IC[nx*N_xdot+nxdot][3] = nat+1;
IC[nx][nxdot] = nat+1;
break;
}
else if(test(yss, Ndt-Nss+1, AT, nat*3, tol)) // if same attractor
{
//IC[nx*N_xdot+nxdot][0] = ic[0]; // fill IC matrix
//IC[nx*N_xdot+nxdot][1] = ic[1];
//IC[nx*N_xdot+nxdot][3] = nat+1;
IC[nx][nxdot] = nat+1;
break;
}
}
}
else // chaotic
{
//IC[nx*N_xdot+nxdot][0] = ic[0]; // fill IC matrix
//IC[nx*N_xdot+nxdot][1] = ic[1];
//IC[nx*N_xdot+nxdot][3] = 0;
IC[nx][nxdot] = 0;
}
for(int i=0; i[] yss[i];}
delete [] yss;
yss = 0;
}
}
delete [] x0;
x0 = 0;
delete [] xdot0;
xdot0 = 0;
for(int i=0; i[] y[i];}
delete [] y;
y = 0;
delete [] Qvolt;
Qvolt = 0;
}
void task_t3(double *d_params, int *i_params, double *param, double** AT, double** IC){
double t0 = d_params[0];
double t1 = d_params[1];
double h = d_params[2];
double tol = d_params[3];
double x0_min = d_params[4];
double x0_max = d_params[5];
double xdot0_min = d_params[6];
double xdot0_max = d_params[7];
double v0 = d_params[8];
int Ndt = i_params[0];
int Nss = i_params[1];
int Nat = i_params[2];
int N_x = i_params[3];
int N_xdot = i_params[4];
// vectors
double *x0 = new double[N_x]; // vector of x0 values
linspace(x0_min, x0_max, N_x, x0);
double *xdot0 = new double[N_xdot]; // vector of xdot0 values
linspace(xdot0_min, xdot0_max, N_xdot, xdot0);
double ic[3];
ic[2] = v0;
// response matrix (dimensions: [Ndt][3])
double **y = new double*[Ndt];
for(int i=0; i[i] = new double[3];}
double *Qvolt = new double[(Ndt-Nss)/50];
// calculating
for(int nx=2; nx[0] = x0[nx]; // updating displacement
for(int nxdot=0; nxdot[1] = xdot0[nxdot]; // updating velocity
RKF4(ic, 3, param, t0, t1, h, pvi, y); // Integrate the dynamical system with RKF4
double **yss = new double *[Ndt-Nss+1]; // steady-state response
for(int i=0; i[i] = new double[3];
yss[i][0] = y[Nss+i-1][0];
yss[i][1] = y[Nss+i-1][1];
yss[i][2] = y[Nss+i-1][2];
}
for(int i=0; i<(Ndt-Nss)/50; i++){
Qvolt[i] = yss[(i*50)][2];
}
if(abs(z1test(Qvolt, (Ndt-Nss)/50)) < 0.2) // if regular
{
for(int nat=0; nat[0][nat*3] == 0) // if attractor is empty
{
fill_AT(yss, 0, Ndt-Nss, AT, nat*3, nat*3+2); // fill attractor
//IC[nx*N_xdot+nxdot][0] = ic[0]; // fill IC matrix
//IC[nx*N_xdot+nxdot][1] = ic[1];
//IC[nx*N_xdot+nxdot][3] = nat+1;
IC[nx][nxdot] = nat+1;
break;
}
else if(test(yss, Ndt-Nss+1, AT, nat*3, tol)) // if same attractor
{
//IC[nx*N_xdot+nxdot][0] = ic[0]; // fill IC matrix
//IC[nx*N_xdot+nxdot][1] = ic[1];
//IC[nx*N_xdot+nxdot][3] = nat+1;
IC[nx][nxdot] = nat+1;
break;
}
}
}
else // chaotic
{
//IC[nx*N_xdot+nxdot][0] = ic[0]; // fill IC matrix
//IC[nx*N_xdot+nxdot][1] = ic[1];
//IC[nx*N_xdot+nxdot][3] = 0;
IC[nx][nxdot] = 0;
}
for(int i=0; i[] yss[i];}
delete [] yss;
yss = 0;
}
}
delete [] x0;
x0 = 0;
delete [] xdot0;
xdot0 = 0;
for(int i=0; i[] y[i];}
delete [] y;
y = 0;
delete [] Qvolt;
Qvolt = 0;
}
void task_t4(double *d_params, int *i_params, double *param, double** AT, double** IC){
double t0 = d_params[0];
double t1 = d_params[1];
double h = d_params[2];
double tol = d_params[3];
double x0_min = d_params[4];
double x0_max = d_params[5];
double xdot0_min = d_params[6];
double xdot0_max = d_params[7];
double v0 = d_params[8];
int Ndt = i_params[0];
int Nss = i_params[1];
int Nat = i_params[2];
int N_x = i_params[3];
int N_xdot = i_params[4];
// vectors
double *x0 = new double[N_x]; // vector of x0 values
linspace(x0_min, x0_max, N_x, x0);
double *xdot0 = new double[N_xdot]; // vector of xdot0 values
linspace(xdot0_min, xdot0_max, N_xdot, xdot0);
double ic[3];
ic[2] = v0;
// response matrix (dimensions: [Ndt][3])
double **y = new double*[Ndt];
for(int i=0; i[i] = new double[3];}
double *Qvolt = new double[(Ndt-Nss)/50];
// calculating
for(int nx=3; nx[0] = x0[nx]; // updating displacement
for(int nxdot=0; nxdot[1] = xdot0[nxdot]; // updating velocity
RKF4(ic, 3, param, t0, t1, h, pvi, y); // Integrate the dynamical system with RKF4
double **yss = new double *[Ndt-Nss+1]; // steady-state response
for(int i=0; i[i] = new double[3];
yss[i][0] = y[Nss+i-1][0];
yss[i][1] = y[Nss+i-1][1];
yss[i][2] = y[Nss+i-1][2];
}
for(int i=0; i<(Ndt-Nss)/50; i++){
Qvolt[i] = yss[(i*50)][2];
}
if(abs(z1test(Qvolt, (Ndt-Nss)/50)) < 0.2) // if regular
{
for(int nat=0; nat[0][nat*3] == 0) // if attractor is empty
{
fill_AT(yss, 0, Ndt-Nss, AT, nat*3, nat*3+2); // fill attractor
//IC[nx*N_xdot+nxdot][0] = ic[0]; // fill IC matrix
//IC[nx*N_xdot+nxdot][1] = ic[1];
//IC[nx*N_xdot+nxdot][3] = nat+1;
IC[nx][nxdot] = nat+1;
break;
}
else if(test(yss, Ndt-Nss+1, AT, nat*3, tol)) // if same attractor
{
//IC[nx*N_xdot+nxdot][0] = ic[0]; // fill IC matrix
//IC[nx*N_xdot+nxdot][1] = ic[1];
//IC[nx*N_xdot+nxdot][3] = nat+1;
IC[nx][nxdot] = nat+1;
break;
}
}
}
else // chaotic
{
//IC[nx*N_xdot+nxdot][0] = ic[0]; // fill IC matrix
//IC[nx*N_xdot+nxdot][1] = ic[1];
//IC[nx*N_xdot+nxdot][3] = 0;
IC[nx][nxdot] = 0;
}
for(int i=0; i[] yss[i];}
delete [] yss;
yss = 0;
}
}
delete [] x0;
x0 = 0;
delete [] xdot0;
xdot0 = 0;
for(int i=0; i[] y[i];}
delete [] y;
y = 0;
delete [] Qvolt;
Qvolt = 0;
}
```c++
// -----------------------------------------------------------------
// header.cpp
// -----------------------------------------------------------------
// This is the function file to compute Basins of attraction and
// others function support for the piezo-magneto-elastic beam
// -----------------------------------------------------------------
// programmer:
// João Peterson (ligier.peterson@gmail.com)
// Americo Cunha (americo.cunha@uerj.br)
//
// last update: Jan, 2021
// -----------------------------------------------------------------
#include // cout
#include // ofstream
#include // round
#include
#include
#include
#include "header.h"
using namespace std;
void percentage(int a, int b){
float n = (float)(a+1) / (float)(b+1);
float percents[10] = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0};
for(int i=0; i<10; i++){
if(n == percents[i]){
cout << "\r";
cout << "[" << (n)*100 << "%]";
}
}
}
void linspace(double a, double b, int N, double output[]){
for(int i=0; i[i] = a + i*(b-a)/(N-1);}
}
bool test(double **yss, int ly, double **AT, int colat, double tol){
bool result = true;
double ATmax, ATmin, ymax, ymin;
for(int i=0; i<3; i++) // changing dim 2
{
ATmax = AT[0][colat+i];
ATmin = AT[0][colat+i];
ymax = yss[0][i];
ymin = yss[0][i];
for(int j=1; j[j][i+colat]);
ATmin = min(ATmin, AT[j][i+colat]);
ymax = max(ymax, yss[j][i]);
ymin = min(ymin, yss[j][i]);
}
if(abs(ATmin-ymin) > tol || abs(ATmax-ymax) > tol)
{
result = false;
}
}
return result;
}
void RKF4(double IC[], int s, double params[], double t0, double t1, double h,
void (*func)(double, double[], double[], double[]), double **y)
{
int Ndt = round((t1-t0)/h); // number of time steps
double *k1, *k2, *k3, *k4; // Runge-Kutta's k
double *y1, *y2, *y3, *y4;
k1 = new double[s];
k2 = new double[s];
k3 = new double[s];
k4 = new double[s];
y1 = new double[s];
y2 = new double[s];
y3 = new double[s];
y4 = new double[s];
double t;
t = t0; // time
for(int i=0; i[0][i] = IC[i];} // initial conditions
for(int i=0; i<(Ndt-1); i++)
{
for(int m=0; m[m] = y[i][m];}
(*func)(t, y1, params, k1);
for(int m=0; m[m] = y[i][m] + 0.5*h*k1[m];}
(*func)(t+0.5*h, y2, params, k2);
for(int m=0; m[m] = y[i][m] + 0.5*h*k2[m];}
(*func)(t+0.5*h, y3, params, k3);
for(int m=0; m[m] = y[i][m] + h*k3[m];}
(*func)(t+h, y4, params, k4);
t += h;
for(int m=0; m[i+1][m] = y[i][m] + h*(0.1666666666666667)*(k1[m] + 2.0*k2[m] + 2.0*k3[m] + k4[m]);
}
}
delete[] k1;
delete[] k2;
delete[] k3;
delete[] k4;
delete[] y1;
delete[] y2;
delete[] y3;
delete[] y4;
}
void pvi(double t, double y[], double param[], double ydot[])
{
// parameters
/*
double ksi = param[0];
double chi = param[1];
double f = param[2];
double omega = param[3];
double lambda = param[4];
double kappa = param[5];
*/
//
// system of equations
ydot[0] = y[1];
ydot[1] = -2.0*param[0]*y[1] + 0.5*y[0]*(1.0 - y[0]*y[0]) + param[1]*y[2] + param[2]*cos(param[3]*t);
//ydot[1] = -2.0*param[0]*y[1] - y[0] + param[1]*y[2] + param[2]*cos(param[3]*t);
ydot[2] = -param[4]*y[2] - param[5]*y[1];
}
// ----------------------------------------------------------------------------
void fill_AT(double **y, int RowYStart, int RowYEnd, double **AT, int ColATStart, int ColATEnd)
{
for(int j=ColATStart; j[i-RowYStart][j] = y[i][j-ColATStart];
}
}
}
// ----------------------------------------------------------------------------
double z1test(double x[], int size_x)
{
int j[size_x];
for(int i=0; i[i] = i+1;}
int roundN10 = round(size_x/10);
double t[roundN10];
for(int i=0; i[i] = i+1;}
double M[roundN10];
int ns = 300; // number of samples
double c[ns];
for(int i=0; i[i] = ((double) rand()/(RAND_MAX))*4*atan(1.0);}
double kcorr[ns];
double *p, *q;
p = new double[size_x];
q = new double[size_x];
double xcos[size_x], xsin[size_x];
for(int its=0; its[i] = x[i]*cos(j[i]*c[its]);
xsin[i] = x[i]*sin(j[i]*c[its]);
}
cumsum(xcos, size_x, p);
cumsum(xsin, size_x, q);
for(int n=0; n[size_x-n+1];
for(int i=0; i<(size_x-n+1); i++){
aux[i] = pow((p[i+n+1]-p[i]), 2.0) + pow((q[i+n+1]-q[i]), 2.0);
}
M[n] = mean(aux, size_x-n+1) - pow(mean(x, size_x), 2.0)*(1.0-cos(n*c[its]))/(1.0-cos(c[its]));
}
kcorr[its] = corr(t, M, roundN10);
}
vector kcorr_less;
vector kcorr_more;
for(int i=0; i[i] < mean(c, ns))
{
kcorr_less.push_back(kcorr[i]);
}
else
{
kcorr_more.push_back(kcorr[i]);
}
}
double array_k_less[kcorr_less.size()];
double array_k_more[kcorr_more.size()];
// checking for oversampling
/*
if( ( (max(x) - min(x))/(mean_abs_diff(x, size_x)) ) > 10 ||
median(array_k_less, int(kcorr_less.size())) - median(array_k_more, int(kcorr_more.size())) > 0.5)
{
cout << "data is probably oversampled." << endl;
cout << "Use coarser sampling or reduce the maximum value of c." << endl;
}
*/
if( ( (*max_element(x, x+size_x) - *min_element(x, x+size_x))/(mean_abs_diff(x, size_x)) ) > 10 ||
median(array_k_less, int(kcorr_less.size())) - median(array_k_more, int(kcorr_more.size())) > 0.5)
{
cout << "data is probably oversampled." << endl;
cout << "Use coarser sampling or reduce the maximum value of c." << endl;
}
delete [] p;
delete [] q;
return median(kcorr, ns);
}
// ----------------------------------------------------------------------------
double mean_abs_diff(double x[], int size_x)
{
double y[size_x-1];
for(int i=0; i[i] = abs(x[i+1] - x[i]);
}
return mean(y, size_x-1);
}
// ----------------------------------------------------------------------------
void cumsum(double x[], int size_x, double result[])
{
result[0] = x[0];
for(int i=1; i[i] = x[i] + result[i-1];
}
}
// ----------------------------------------------------------------------------
double sum(double a[], int size_a)
{
double s = 0;
for (int i = 0; i < size_a; i++)
{
s += a[i];
}
return s;
}
// ----------------------------------------------------------------------------
double mean(double a[], int size_a)
{
return sum(a, size_a) / size_a;
}
// ----------------------------------------------------------------------------
double corr(double x[], double y[], int size_n)
{
/* pearson coefficient
_n
\ (x_i - mean(x))(y_i - mean(y))
/_1
r = ------------------------------------------------------------
sqrt(sum((x_i - mean(x))^2)) * sqrt(sum((y_i - mean(y))^2))
*/
double N = 0;
double D1 = 0;
double D2 = 0;
double mean_x = mean(x, size_n);
double mean_y = mean(y, size_n);
for(int i=0; i[i] - mean_x)*(y[i] - mean_y);
D1 += pow((x[i] - mean_x), 2.0);
D2 += pow((y[i] - mean_y), 2.0);
}
D1 = sqrt(D1);
D2 = sqrt(D2);
return (N/(D1*D2));
}
// ----------------------------------------------------------------------------
double median(double x[], int size_x)
{
double result;
my_sort(x, size_x);
if(size_x % 2 == 0)
result = (x[(size_x-1)/2] + x[(size_x+1)/2])/2;
else
result = x[size_x/2];
return result;
}
// ----------------------------------------------------------------------------
void my_sort(double x[], int size_x)
{
double a;
for(int i=0; i[j][i]){
a = x[j];
x[j] = x[i];
x[i] = a;
}
}
}
}
c++
/* -----------------------------------------------------------------
* header.h
* -----------------------------------------------------------------
* This is the function file to compute Basins of attraction and
* others function support for the piezo-magneto-elastic beam
* -----------------------------------------------------------------
* programmer:
* Jo?o Peterson (ligier.peterson@gmail.com)
* Americo Cunha (americo.cunha@uerj.br)
*
* last update: Jan, 2021
*/ -----------------------------------------------------------------
#ifndef HEADER_H_INCLUDED
#define HEADER_H_INCLUDED
#include // cout
#include // ofstream
#include // round
#include
#include
using namespace std;
void percentage(int a, int b);
void linspace(double a, double b, int N, double output[]);
bool test(double **yss, int ly, double **AT, int colat, double tol);
void RKF4(double IC[], int s, double params[], double t0, double t1, double h,
void (*func)(double, double[], double[], double[]), double **y);
void pvi(double t, double y[], double param[], double ydot[]);
void fill_AT(double **y, int RowYStart, int RowYEnd, double **AT, int ColATStart, int ColATEnd);
double z1test(double x[], int size_x);
double mean_abs_diff(double x[], int size_x);
void cumsum(double x[], int size_x, double result[]);
double sum(double a[], int size_a);
double mean(double a[], int size_a);
double corr(double x[], double y[], int size_n);
double median(double x[], int size_x);
void my_sort(double x[], int size_x);
#endif // HEADER_H_INCLUDED