add old lab3 and basic szamszim
This commit is contained in:
parent
151fb63720
commit
12e8482055
|
|
@ -0,0 +1,34 @@
|
|||
#+OPTIONS: tex:t
|
||||
#+AUTHOR: Rédl Anna, Barna Zsombor
|
||||
#+DATE: 2023. 07. 25.
|
||||
#+TITLE: Összegyűjtött
|
||||
#+LATEX_CLASS: article
|
||||
#+LATEX_CLASS_OPTIONS: [12pt,a4paper]
|
||||
#+LATEX_HEADER: \usepackage[utf8]{inputenc}
|
||||
#+LATEX_HEADER: \usepackage[magyar]{babel}
|
||||
#+LATEX_HEADER: \hypersetup{ colorlinks=true,linkcolor=blue, linktoc=all, filecolor=magenta, urlcolor=cyan, pdfstartview=FitB,}
|
||||
#+LATEX_HEADER: \usepackage{multirow}
|
||||
#+LATEX_HEADER: \usepackage{braket}
|
||||
#+LATEX_HEADER: \urlstyle{same}
|
||||
#+OPTIONS: toc:nil
|
||||
|
||||
#+tblname: deuterium
|
||||
| $\nu_D[cm^{-1}]$ | $I_D[km/mol]$ |
|
||||
|------------------+---------------|
|
||||
| 361.53 | 0.000 |
|
||||
| 362.23 | 0.000 |
|
||||
| 592.82 | 0.000 |
|
||||
| 592.85 | 0.000 |
|
||||
| 615.29 | 0.000 |
|
||||
#+call: plot(inputs=deuterium, fname="fig2.png") :noexport:
|
||||
|
||||
|
||||
#+name: plot
|
||||
#+begin_src python :var inputs=benzene :var fname="" :noweb yes
|
||||
func_list = []
|
||||
for row in inputs:
|
||||
nu = row[0]
|
||||
A = row[1]
|
||||
if A > 0.001:
|
||||
func_list.append(get_distr(A, nu))
|
||||
#+end_src
|
||||
|
|
@ -0,0 +1,36 @@
|
|||
##################################################
|
||||
# Author: József Stéger
|
||||
# Date: 02 13 2017
|
||||
# Usage (tested in linux):
|
||||
# make
|
||||
# Summary:
|
||||
# 1. compiles objects in ./build/
|
||||
# 2. creates lib archive in .
|
||||
##################################################
|
||||
|
||||
UTILS=datafit fft interpolate matrix odeint optimize random vector
|
||||
BUILD=build
|
||||
|
||||
OBJ_LIBS=$(patsubst %,$(BUILD)/%.o,$(UTILS))
|
||||
A_LIB=libcpl.a
|
||||
|
||||
.PHONY: all distclean clean
|
||||
|
||||
all: $(A_LIB)
|
||||
|
||||
clean:
|
||||
rm -rf $(A_LIB) $(BUILD)
|
||||
|
||||
distclean: clean
|
||||
rm -rf $(PDIR)
|
||||
|
||||
$(OBJ_LIBS):
|
||||
@make $(BUILD)
|
||||
g++ -Wall -c -I. -O2 $(patsubst $(BUILD)/%.o,./%.cpp,$@) -o $@
|
||||
|
||||
$(A_LIB): $(OBJ_LIBS)
|
||||
ar -srcl $@ $?
|
||||
|
||||
$(BUILD):
|
||||
@mkdir -p $@
|
||||
|
||||
Binary file not shown.
|
|
@ -0,0 +1,178 @@
|
|||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <cstdlib>
|
||||
#include <iostream>
|
||||
#include <limits>
|
||||
#include <vector>
|
||||
using namespace std;
|
||||
|
||||
#include "vector.hpp"
|
||||
#include "datafit.hpp"
|
||||
|
||||
namespace cpl {
|
||||
|
||||
inline double SQR(double a) {
|
||||
return a * a;
|
||||
}
|
||||
|
||||
static void error(const char *msg) {
|
||||
cerr << msg << endl;
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
double gammln(const double xx)
|
||||
{
|
||||
int j;
|
||||
double x,y,tmp,ser;
|
||||
static const double cof[6]=
|
||||
{76.18009172947146,-86.50532032941677,
|
||||
24.01409824083091,-1.231739572450155,0.1208650973866179e-2,
|
||||
-0.5395239384953e-5};
|
||||
|
||||
y=x=xx;
|
||||
tmp=x+5.5;
|
||||
tmp -= (x+0.5)*log(tmp);
|
||||
ser=1.000000000190015;
|
||||
for (j=0;j<6;j++) ser += cof[j]/++y;
|
||||
return -tmp+log(2.5066282746310005*ser/x);
|
||||
}
|
||||
|
||||
void gser(double& gamser, const double a, const double x, double& gln)
|
||||
{
|
||||
const int ITMAX=100;
|
||||
const double EPS=numeric_limits<double>::epsilon();
|
||||
int n;
|
||||
double sum,del,ap;
|
||||
|
||||
gln=gammln(a);
|
||||
if (x <= 0.0) {
|
||||
if (x < 0.0) error("x less than 0 in routine gser");
|
||||
gamser=0.0;
|
||||
return;
|
||||
} else {
|
||||
ap=a;
|
||||
del=sum=1.0/a;
|
||||
for (n=0;n<ITMAX;n++) {
|
||||
++ap;
|
||||
del *= x/ap;
|
||||
sum += del;
|
||||
if (abs(del) < abs(sum)*EPS) {
|
||||
gamser=sum*exp(-x+a*log(x)-gln);
|
||||
return;
|
||||
}
|
||||
}
|
||||
error("a too large, ITMAX too small in routine gser");
|
||||
return;
|
||||
}
|
||||
}
|
||||
|
||||
void gcf(double& gammcf, const double a, const double x, double& gln)
|
||||
{
|
||||
const int ITMAX=100;
|
||||
const double EPS=numeric_limits<double>::epsilon();
|
||||
const double FPMIN=numeric_limits<double>::min()/EPS;
|
||||
int i;
|
||||
double an,b,c,d,del,h;
|
||||
|
||||
gln=gammln(a);
|
||||
b=x+1.0-a;
|
||||
c=1.0/FPMIN;
|
||||
d=1.0/b;
|
||||
h=d;
|
||||
for (i=1;i<=ITMAX;i++) {
|
||||
an = -i*(i-a);
|
||||
b += 2.0;
|
||||
d=an*d+b;
|
||||
if (abs(d) < FPMIN) d=FPMIN;
|
||||
c=b+an/c;
|
||||
if (abs(c) < FPMIN) c=FPMIN;
|
||||
d=1.0/d;
|
||||
del=d*c;
|
||||
h *= del;
|
||||
if (abs(del-1.0) <= EPS) break;
|
||||
}
|
||||
if (i > ITMAX) error("a too large, ITMAX too small in gcf");
|
||||
gammcf=exp(-x+a*log(x)-gln)*h;
|
||||
}
|
||||
|
||||
double gammq(const double a, const double x)
|
||||
{
|
||||
double gamser,gammcf,gln;
|
||||
|
||||
if (x < 0.0 || a <= 0.0)
|
||||
error("Invalid arguments in routine gammq");
|
||||
if (x < a+1.0) {
|
||||
gser(gamser,a,x,gln);
|
||||
return 1.0-gamser;
|
||||
} else {
|
||||
gcf(gammcf,a,x,gln);
|
||||
return gammcf;
|
||||
}
|
||||
}
|
||||
|
||||
void fit(
|
||||
const Vector& x,
|
||||
const Vector& y,
|
||||
const Vector& sig,
|
||||
const bool mwt,
|
||||
double& a,
|
||||
double& b,
|
||||
double& siga,
|
||||
double& sigb,
|
||||
double& chi2,
|
||||
double& q)
|
||||
{
|
||||
int i;
|
||||
double wt,t,sxoss,sx=0.0,sy=0.0,st2=0.0,ss,sigdat;
|
||||
|
||||
int ndata=x.size();
|
||||
b=0.0;
|
||||
if (mwt) {
|
||||
ss=0.0;
|
||||
for (i=0;i<ndata;i++) {
|
||||
wt=1.0/SQR(sig[i]);
|
||||
ss += wt;
|
||||
sx += x[i]*wt;
|
||||
sy += y[i]*wt;
|
||||
}
|
||||
} else {
|
||||
for (i=0;i<ndata;i++) {
|
||||
sx += x[i];
|
||||
sy += y[i];
|
||||
}
|
||||
ss=ndata;
|
||||
}
|
||||
sxoss=sx/ss;
|
||||
if (mwt) {
|
||||
for (i=0;i<ndata;i++) {
|
||||
t=(x[i]-sxoss)/sig[i];
|
||||
st2 += t*t;
|
||||
b += t*y[i]/sig[i];
|
||||
}
|
||||
} else {
|
||||
for (i=0;i<ndata;i++) {
|
||||
t=x[i]-sxoss;
|
||||
st2 += t*t;
|
||||
b += t*y[i];
|
||||
}
|
||||
}
|
||||
b /= st2;
|
||||
a=(sy-sx*b)/ss;
|
||||
siga=sqrt((1.0+sx*sx/(ss*st2))/ss);
|
||||
sigb=sqrt(1.0/st2);
|
||||
chi2=0.0;
|
||||
q=1.0;
|
||||
if (!mwt) {
|
||||
for (i=0;i<ndata;i++)
|
||||
chi2 += SQR(y[i]-a-b*x[i]);
|
||||
sigdat=sqrt(chi2/(ndata-2));
|
||||
siga *= sigdat;
|
||||
sigb *= sigdat;
|
||||
} else {
|
||||
for (i=0;i<ndata;i++)
|
||||
chi2 += SQR((y[i]-a-b*x[i])/sig[i]);
|
||||
if (ndata>2) q=gammq(0.5*(ndata-2),0.5*chi2);
|
||||
}
|
||||
}
|
||||
|
||||
} /* end namespace cpl */
|
||||
|
|
@ -0,0 +1,23 @@
|
|||
#ifndef CPL_DATAFIT_HPP
|
||||
#define CPL_DATAFIT_HPP
|
||||
|
||||
namespace cpl {
|
||||
|
||||
class Vector;
|
||||
|
||||
extern void fit( // Least Squares fit to y(x) = a + b x
|
||||
const Vector& x, // Input: Vector of x values
|
||||
const Vector& y, // Input: Vector of y values
|
||||
const Vector& sig, // Input: Vector of individual standard deviations
|
||||
const bool mwt, // Input: if false sig's assumed to be unavailable
|
||||
double& a, // Output: intercept a
|
||||
double& b, // Output: slope b
|
||||
double& siga, // Output: uncertainty in a
|
||||
double& sigb, // Output: uncertainty in b
|
||||
double& chi2, // Output: chi-square
|
||||
double& q // Output: goodness of fit
|
||||
);
|
||||
|
||||
} /* end namespace cpl */
|
||||
|
||||
#endif /* CPL_DATAFIT_HPP */
|
||||
|
|
@ -0,0 +1,214 @@
|
|||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <cstdlib>
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
using namespace std;
|
||||
|
||||
#include "fft.hpp"
|
||||
#include "vector.hpp"
|
||||
|
||||
namespace cpl {
|
||||
|
||||
ComplexVector::ComplexVector(int dim) {
|
||||
v = new std::complex<double> [this->dim = dim];
|
||||
for (int i = 0; i < dim; i++) v[i] = 0.0;
|
||||
}
|
||||
|
||||
ComplexVector::ComplexVector(const ComplexVector& cv) {
|
||||
v = new std::complex<double> [dim = cv.dim];
|
||||
for (int i = 0; i < dim; i++) v[i] = cv.v[i];
|
||||
}
|
||||
|
||||
ComplexVector& ComplexVector::operator = (const ComplexVector& cv) {
|
||||
if (this != &cv) {
|
||||
if (dim != cv.dim) {
|
||||
delete [] v;
|
||||
v = new std::complex<double> [dim = cv.dim];
|
||||
}
|
||||
for (int i = 0; i < dim; i++) v[i] = cv[i];
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
// FFT implementation
|
||||
|
||||
void FFT::transform(ComplexVector& data) {
|
||||
N = data.dimension();
|
||||
f = &data;
|
||||
bitReverse();
|
||||
for (int n = 1; n < N; n *= 2)
|
||||
DanielsonLanczos(n);
|
||||
for (int i = 0; i < N; ++i)
|
||||
(*f)[i] /= std::sqrt(double(N));
|
||||
}
|
||||
|
||||
void FFT::inverseTransform(ComplexVector& data) {
|
||||
inverse = true;
|
||||
transform(data);
|
||||
inverse = false;
|
||||
}
|
||||
|
||||
void FFT::bitReverse() {
|
||||
int j = 1;
|
||||
for (int i = 1; i < N; ++i) {
|
||||
if (i < j) {
|
||||
std::complex<double> temp = (*f)[i-1];
|
||||
(*f)[i-1] = (*f)[j-1];
|
||||
(*f)[j-1] = temp;
|
||||
}
|
||||
int k = N / 2;
|
||||
while ( k < j ) {
|
||||
j -= k;
|
||||
k /= 2;
|
||||
}
|
||||
j += k;
|
||||
}
|
||||
}
|
||||
|
||||
void FFT::DanielsonLanczos(int n) {
|
||||
const double pi = 4 * atan(1.0);
|
||||
std::complex<double> W(0, pi / n);
|
||||
W = inverse ? std::exp(W) : std::exp(-W);
|
||||
std::complex<double> W_j(1, 0);
|
||||
for (int j = 0; j < n; ++j) {
|
||||
for (int i = j; i < N; i += 2 * n) {
|
||||
std::complex<double> temp = W_j * (*f)[n+i];
|
||||
(*f)[n+i] = (*f)[i] - temp;
|
||||
(*f)[i] += temp;
|
||||
}
|
||||
W_j *= W;
|
||||
}
|
||||
}
|
||||
|
||||
Vector FFT::power(ComplexVector& data) {
|
||||
Vector P(1 + N / 2);
|
||||
P[0] = std::norm(data[0]) / double(N);
|
||||
for (int i = 1; i < N / 2; i++)
|
||||
P[i] = (std::norm(data[i]) + std::norm(data[N-i])) / double(N);
|
||||
P[N/2] = std::norm(data[N/2]) / double(N);
|
||||
return P;
|
||||
}
|
||||
|
||||
static void error(const std::string str) {
|
||||
std::cerr << "cpl::ComplexMatrix error: " << str << std::endl;
|
||||
std::exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
static void error(const std::string str, int i) {
|
||||
std::ostringstream os;
|
||||
os << str << " " << i;
|
||||
error(os.str());
|
||||
}
|
||||
|
||||
ComplexMatrix::ComplexMatrix(int rows, int cols, complex<double> d) {
|
||||
if (rows <= 0)
|
||||
error("bad number of rows", rows);
|
||||
if (cols <= 0)
|
||||
error("bad number of columns", cols);
|
||||
this->rows = rows;
|
||||
this->cols = cols;
|
||||
m = new complex<double>* [rows];
|
||||
for (int i = 0; i < rows; i++) {
|
||||
m[i] = new complex<double> [cols];
|
||||
for (int j = 0; j < cols; j++)
|
||||
m[i][j] = d;
|
||||
}
|
||||
}
|
||||
|
||||
ComplexMatrix::ComplexMatrix(const ComplexMatrix& mat) {
|
||||
rows = mat.rows;
|
||||
cols = mat.cols;
|
||||
|
||||
m = new complex<double>* [rows];
|
||||
for (int i = 0; i < rows; i++) {
|
||||
m[i] = new complex<double> [cols];
|
||||
for (int j = 0; j < cols; j++)
|
||||
m[i][j] = mat[i][j];
|
||||
}
|
||||
}
|
||||
|
||||
ComplexMatrix::~ComplexMatrix() {
|
||||
for (int i = 0; i < rows; i++)
|
||||
delete [] m[i];
|
||||
delete [] m;
|
||||
}
|
||||
|
||||
const complex<double>* ComplexMatrix::operator[](const int row) const {
|
||||
if (row < 0 || row >= rows)
|
||||
error("bad row index", row);
|
||||
return m[row];
|
||||
}
|
||||
|
||||
complex<double>* ComplexMatrix::operator[](const int row) {
|
||||
if (row < 0 || row >= rows)
|
||||
error("bad row index", row);
|
||||
return m[row];
|
||||
}
|
||||
|
||||
ComplexMatrix& ComplexMatrix::operator=(const ComplexMatrix& mat) {
|
||||
if (this != &mat) {
|
||||
if (rows != mat.rows || cols != mat.cols) {
|
||||
for (int i = 0; i < rows; i++)
|
||||
delete [] m[i];
|
||||
delete [] m;
|
||||
rows = mat.rows;
|
||||
cols = mat.cols;
|
||||
m = new complex<double>* [rows];
|
||||
for (int i = 0; i < rows; i++)
|
||||
m[i] = new complex<double> [cols];
|
||||
}
|
||||
for (int i = 0; i < rows; i++)
|
||||
for (int j = 0; j < cols; j++)
|
||||
m[i][j] = mat[i][j];
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
void FFT2::transform(ComplexMatrix& data) {
|
||||
FFT fft;
|
||||
int rows = data.numRows();
|
||||
int cols = data.numCols();
|
||||
ComplexVector r(cols), c(rows);
|
||||
// FFT rows of data
|
||||
for (int j = 0; j < rows; j++) {
|
||||
for (int k = 0; k < cols; k++)
|
||||
r[k] = data[j][k];
|
||||
fft.transform(r);
|
||||
for (int k = 0; k < cols; k++)
|
||||
data[j][k] = r[k];
|
||||
}
|
||||
// FFT columns of data
|
||||
for (int k = 0; k < cols; k++) {
|
||||
for (int j = 0; j < rows; j++)
|
||||
c[j] = data[j][k];
|
||||
fft.transform(c);
|
||||
for (int j = 0; j < rows; j++)
|
||||
data[j][k] = c[j];
|
||||
}
|
||||
}
|
||||
|
||||
void FFT2::inverseTransform(ComplexMatrix& data) {
|
||||
FFT fft;
|
||||
int rows = data.numRows();
|
||||
int cols = data.numCols();
|
||||
ComplexVector r(cols), c(rows);
|
||||
// FFT rows of data
|
||||
for (int j = 0; j < rows; j++) {
|
||||
for (int k = 0; k < cols; k++)
|
||||
r[k] = data[j][k];
|
||||
fft.inverseTransform(r);
|
||||
for (int k = 0; k < cols; k++)
|
||||
data[j][k] = r[k];
|
||||
}
|
||||
// FFT columns of data
|
||||
for (int k = 0; k < cols; k++) {
|
||||
for (int j = 0; j < rows; j++)
|
||||
c[j] = data[j][k];
|
||||
fft.inverseTransform(c);
|
||||
for (int j = 0; j < rows; j++)
|
||||
data[j][k] = c[j];
|
||||
}
|
||||
}
|
||||
|
||||
} /* end namespace cpl */
|
||||
|
|
@ -0,0 +1,75 @@
|
|||
#ifndef CPL_FFT_HPP
|
||||
#define CPL_FFT_HPP
|
||||
|
||||
#include <complex>
|
||||
#include <iostream>
|
||||
|
||||
namespace cpl {
|
||||
|
||||
class Vector;
|
||||
|
||||
class ComplexVector {
|
||||
public:
|
||||
ComplexVector(int dim = 1);
|
||||
ComplexVector(const ComplexVector& cv);
|
||||
~ComplexVector() { delete [] v; }
|
||||
|
||||
int dimension() const { return dim; }
|
||||
const std::complex<double> operator[](const int i) const { return v[i]; }
|
||||
std::complex<double>& operator[](const int i) { return v[i]; }
|
||||
ComplexVector& operator = (const ComplexVector& cv);
|
||||
|
||||
private:
|
||||
int dim;
|
||||
std::complex<double> *v;
|
||||
};
|
||||
|
||||
class FFT {
|
||||
public:
|
||||
FFT() { N = 0; f = 0; inverse = false; }
|
||||
void transform(ComplexVector& data);
|
||||
void inverseTransform(ComplexVector& data);
|
||||
Vector power(ComplexVector& data);
|
||||
|
||||
private:
|
||||
int N;
|
||||
ComplexVector *f;
|
||||
bool inverse;
|
||||
|
||||
void bitReverse();
|
||||
void DanielsonLanczos(int n);
|
||||
};
|
||||
|
||||
class ComplexMatrix {
|
||||
public:
|
||||
|
||||
ComplexMatrix(int rows=1, int cols=1, complex<double> d=0);
|
||||
ComplexMatrix(const ComplexMatrix& m);
|
||||
~ComplexMatrix();
|
||||
|
||||
int numRows() const { return rows; }
|
||||
int numCols() const { return cols; }
|
||||
const std::complex<double>* operator[](const int row) const;
|
||||
std::complex<double>* operator[](const int row);
|
||||
std::complex<double>& operator()(int row, int col);
|
||||
ComplexMatrix& operator=(const ComplexMatrix& m);
|
||||
|
||||
private:
|
||||
int rows;
|
||||
int cols;
|
||||
std::complex<double> **m;
|
||||
};
|
||||
|
||||
class FFT2 {
|
||||
public:
|
||||
FFT2() { }
|
||||
void transform(ComplexMatrix& data);
|
||||
void inverseTransform(ComplexMatrix& data);
|
||||
|
||||
private:
|
||||
|
||||
};
|
||||
|
||||
} /* end namespace cpl */
|
||||
|
||||
#endif /* CPL_FFT_HPP */
|
||||
|
|
@ -0,0 +1,59 @@
|
|||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <cstdlib>
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
using namespace std;
|
||||
|
||||
#include "interpolate.hpp"
|
||||
|
||||
namespace cpl {
|
||||
|
||||
// polynomial interpolation routine adapted from Numerical Recipes
|
||||
// polynomial degree inferred from dimension of input vector
|
||||
|
||||
void Interpolator::polint( // uses Neville's algorithm
|
||||
const vector<double>& xa, // vector of abscissa values
|
||||
const vector<double>& ya, // vector of ordinate values
|
||||
const double x, // point x
|
||||
double& y, // interpolated y(x)
|
||||
double& dy) // estimate of error in interpolation
|
||||
{
|
||||
int i,m,ns=0;
|
||||
double den,dif,dift,ho,hp,w;
|
||||
|
||||
int n=xa.size();
|
||||
vector<double> c(n),d(n);
|
||||
dif=abs(x-xa[0]);
|
||||
for (i=0;i<n;i++) {
|
||||
if ((dift=abs(x-xa[i])) < dif) {
|
||||
ns=i;
|
||||
dif=dift;
|
||||
}
|
||||
c[i]=ya[i];
|
||||
d[i]=ya[i];
|
||||
}
|
||||
y=ya[ns--];
|
||||
for (m=1;m<n;m++) {
|
||||
for (i=0;i<n-m;i++) {
|
||||
ho=xa[i]-x;
|
||||
hp=xa[i+m]-x;
|
||||
w=c[i+1]-d[i];
|
||||
if ((den=ho-hp) == 0.0) { // nrerror("Error in routine polint");
|
||||
cerr << "Error in Interpolator::polint\n"
|
||||
<< "Neville's algorithm breaks down when two "
|
||||
<< "input x values are equal within roundoff\n"
|
||||
<< "x, y values" << endl;
|
||||
for (int j = 0; j < n; j++)
|
||||
cerr << xa[j] << '\t' << ya[j] << '\n';
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
den=w/den;
|
||||
d[i]=hp*den;
|
||||
c[i]=ho*den;
|
||||
}
|
||||
y += (dy=(2*(ns+1) < (n-m) ? c[ns+1] : d[ns--]));
|
||||
}
|
||||
}
|
||||
|
||||
} /* end namespace cpl */
|
||||
|
|
@ -0,0 +1,19 @@
|
|||
#ifndef CPL_INTERPOLATE_HPP
|
||||
#define CPL_INTERPOLATE_HPP
|
||||
|
||||
namespace cpl {
|
||||
|
||||
class Interpolator {
|
||||
public:
|
||||
|
||||
static void polint( // adapted from Numerical Recipes
|
||||
const vector<double>& xa, // vector of x values - input
|
||||
const vector<double>& ya, // vector of y values - input
|
||||
const double x, // x at which to find y - input
|
||||
double& y, // value of y - output
|
||||
double& dy); // estimate or error in y
|
||||
};
|
||||
|
||||
} /* end namespace cpl */
|
||||
|
||||
#endif /* CPL_INTERPOLATE_HPP */
|
||||
|
|
@ -0,0 +1,895 @@
|
|||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <cstdlib>
|
||||
#include <iostream>
|
||||
#include <limits>
|
||||
#include <string>
|
||||
#include <sstream>
|
||||
|
||||
#include "matrix.hpp"
|
||||
#include "vector.hpp"
|
||||
|
||||
namespace cpl {
|
||||
|
||||
static void error(const std::string str) {
|
||||
std::cerr << "cpl::Matrix error: " << str << std::endl;
|
||||
std::exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
static void error(const std::string str, int i) {
|
||||
std::ostringstream os;
|
||||
os << str << " " << i;
|
||||
error(os.str());
|
||||
}
|
||||
|
||||
Matrix::Matrix(int rows, int cols, double d) {
|
||||
if (rows <= 0)
|
||||
error("bad number of rows", rows);
|
||||
if (cols <= 0)
|
||||
error("bad number of columns", cols);
|
||||
this->rows = rows;
|
||||
this->cols = cols;
|
||||
m = new double* [rows];
|
||||
for (int i = 0; i < rows; i++) {
|
||||
m[i] = new double [cols];
|
||||
for (int j = 0; j < cols; j++)
|
||||
m[i][j] = d;
|
||||
}
|
||||
}
|
||||
|
||||
Matrix::Matrix(const Matrix& mat) {
|
||||
rows = mat.rows;
|
||||
cols = mat.cols;
|
||||
|
||||
m = new double* [rows];
|
||||
for (int i = 0; i < rows; i++) {
|
||||
m[i] = new double [cols];
|
||||
for (int j = 0; j < cols; j++)
|
||||
m[i][j] = mat[i][j];
|
||||
}
|
||||
}
|
||||
|
||||
Matrix::~Matrix() {
|
||||
for (int i = 0; i < rows; i++)
|
||||
delete [] m[i];
|
||||
delete [] m;
|
||||
}
|
||||
|
||||
const double* Matrix::operator[](const int row) const {
|
||||
if (row < 0 || row >= rows)
|
||||
error("bad row index", row);
|
||||
return m[row];
|
||||
}
|
||||
|
||||
double* Matrix::operator[](const int row) {
|
||||
if (row < 0 || row >= rows)
|
||||
error("bad row index", row);
|
||||
return m[row];
|
||||
}
|
||||
|
||||
Matrix& Matrix::operator=(const Matrix& mat) {
|
||||
if (this != &mat) {
|
||||
if (rows != mat.rows || cols != mat.cols) {
|
||||
for (int i = 0; i < rows; i++)
|
||||
delete [] m[i];
|
||||
delete [] m;
|
||||
rows = mat.rows;
|
||||
cols = mat.cols;
|
||||
m = new double* [rows];
|
||||
for (int i = 0; i < rows; i++)
|
||||
m[i] = new double [cols];
|
||||
}
|
||||
for (int i = 0; i < rows; i++)
|
||||
for (int j = 0; j < cols; j++)
|
||||
m[i][j] = mat[i][j];
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
Matrix& Matrix::operator=(const double d) {
|
||||
for (int i = 0; i < rows; i++)
|
||||
for (int j = 0; j < cols; j++)
|
||||
m[i][j] = d;
|
||||
return *this;
|
||||
}
|
||||
|
||||
Matrix& Matrix::operator+=(const Matrix& mat) {
|
||||
if (this != &mat) {
|
||||
if (rows != mat.rows || cols != mat.cols)
|
||||
error("matrix dimension mismatch");
|
||||
for (int i = 0; i < rows; i++)
|
||||
for (int j = 0; j < cols; j++)
|
||||
m[i][j] += mat[i][j];
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
Matrix& Matrix::operator-=(const Matrix& mat) {
|
||||
if (this != &mat) {
|
||||
if (rows != mat.rows || cols != mat.cols)
|
||||
error("matrix dimension mismatch");
|
||||
for (int i = 0; i < rows; i++)
|
||||
for (int j = 0; j < cols; j++)
|
||||
m[i][j] -= mat[i][j];
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
Matrix& Matrix::operator*=(const double d) {
|
||||
for (int i = 0; i < rows; i++)
|
||||
for (int j = 0; j < cols; j++)
|
||||
m[i][j] *= d;
|
||||
return *this;
|
||||
}
|
||||
|
||||
Matrix& Matrix::operator/=(const double d) {
|
||||
for (int i = 0; i < rows; i++)
|
||||
for (int j = 0; j < cols; j++)
|
||||
m[i][j] /= d;
|
||||
return *this;
|
||||
}
|
||||
|
||||
Matrix operator - (const Matrix& mat) {
|
||||
int rows = mat.numRows();
|
||||
int cols = mat.numCols();
|
||||
Matrix temp(rows, cols);
|
||||
for (int i = 0; i < rows; i++)
|
||||
for (int j = 0; j < cols; j++)
|
||||
temp[i][j] = -mat[i][j];
|
||||
return temp;
|
||||
}
|
||||
|
||||
Matrix operator * (const Matrix& mat, double d) {
|
||||
int rows = mat.numRows();
|
||||
int cols = mat.numCols();
|
||||
Matrix temp(rows, cols);
|
||||
for (int i = 0; i < rows; i++)
|
||||
for (int j = 0; j < cols; j++)
|
||||
temp[i][j] = mat[i][j] * d;
|
||||
return temp;
|
||||
}
|
||||
|
||||
Matrix operator * (double d, const Matrix& mat) {
|
||||
int rows = mat.numRows();
|
||||
int cols = mat.numCols();
|
||||
Matrix temp(rows, cols);
|
||||
for (int i = 0; i < rows; i++)
|
||||
for (int j = 0; j < cols; j++)
|
||||
temp[i][j] = d * mat[i][j];
|
||||
return temp;
|
||||
}
|
||||
|
||||
Matrix operator / (const Matrix& mat, double d) {
|
||||
int rows = mat.numRows();
|
||||
int cols = mat.numCols();
|
||||
Matrix temp(rows, cols);
|
||||
for (int i = 0; i < rows; i++)
|
||||
for (int j = 0; j < cols; j++)
|
||||
temp[i][j] = mat[i][j] / d;
|
||||
return temp;
|
||||
}
|
||||
|
||||
Matrix operator + (const Matrix& m1, const Matrix& m2) {
|
||||
int rows = m1.numRows();
|
||||
int cols = m1.numCols();
|
||||
if (rows != m2.numRows() || cols != m2.numCols())
|
||||
error("matrix dimension mismatch");
|
||||
Matrix temp(rows, cols);
|
||||
for (int i = 0; i < rows; i++)
|
||||
for (int j = 0; j < cols; j++)
|
||||
temp[i][j] = m1[i][j] + m2[i][j];
|
||||
return temp;
|
||||
}
|
||||
|
||||
Matrix operator - (const Matrix& m1, const Matrix& m2) {
|
||||
int rows = m1.numRows();
|
||||
int cols = m1.numCols();
|
||||
if (rows != m2.numRows() || cols != m2.numCols())
|
||||
error("matrix dimension mismatch");
|
||||
Matrix temp(rows, cols);
|
||||
for (int i = 0; i < rows; i++)
|
||||
for (int j = 0; j < cols; j++)
|
||||
temp[i][j] = m1[i][j] - m2[i][j];
|
||||
return temp;
|
||||
}
|
||||
|
||||
Matrix operator * (const Matrix& m1, const Matrix& m2) {
|
||||
int n = m1.numCols();
|
||||
if (n != m2.numRows())
|
||||
error("matrix dimension mismatch");
|
||||
int rows = m1.numRows();
|
||||
int cols = m2.numCols();
|
||||
Matrix temp(rows, cols);
|
||||
for (int i = 0; i < rows; i++)
|
||||
for (int j = 0; j < cols; j++)
|
||||
for (int k = 0; k < n; k++)
|
||||
temp[i][j] += m1[i][k] * m2[k][j];
|
||||
return temp;
|
||||
}
|
||||
|
||||
Matrix Matrix::transpose() {
|
||||
Matrix temp(cols, rows);
|
||||
for (int i = 0; i < rows; i++)
|
||||
for (int j = 0; j < cols; j++)
|
||||
temp[j][i] = m[i][j];
|
||||
return temp;
|
||||
}
|
||||
|
||||
std::ostream& operator<<(std::ostream& os, const Matrix& mat) {
|
||||
for (int i = 0; i < mat.rows; i++) {
|
||||
for (int j = 0; j < mat.cols; j++) {
|
||||
os << '\t' << mat.m[i][j];
|
||||
}
|
||||
os << '\n';
|
||||
}
|
||||
return os;
|
||||
}
|
||||
|
||||
void solveGaussJordan(Matrix& A, Matrix& B) {
|
||||
|
||||
int n = A.numRows();
|
||||
if (A.numCols() != n)
|
||||
error("matrix dimension mismatch");
|
||||
if (B.numRows() != n)
|
||||
error("matrix dimension mismatch");
|
||||
int m = B.numCols();
|
||||
|
||||
int *indxc = new int [n];
|
||||
int *indxr = new int [n];
|
||||
int *ipiv = new int [n];
|
||||
|
||||
for (int j = 0; j < n; j++)
|
||||
ipiv[j] = 0;
|
||||
|
||||
for (int i = 0; i < n; i++) {
|
||||
double big = 0;
|
||||
int irow, icol;
|
||||
for (int j = 0; j < n; j++) {
|
||||
if (ipiv[j] != 1) {
|
||||
for (int k = 0; k < n; k++) {
|
||||
if (ipiv[k] == 0) {
|
||||
double Ajk = A[j][k];
|
||||
if (std::abs(Ajk) >= big) {
|
||||
big = std::abs(Ajk);
|
||||
irow = j;
|
||||
icol = k;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
++ipiv[icol];
|
||||
|
||||
if (irow != icol) {
|
||||
for (int j = 0; j < n; j++)
|
||||
std::swap(A[irow][j], A[icol][j]);
|
||||
for (int j = 0; j < m; j++)
|
||||
std::swap(B[irow][j], B[icol][j]);
|
||||
}
|
||||
|
||||
indxr[i] = irow;
|
||||
indxc[i] = icol;
|
||||
if (A[icol][icol] == 0.0)
|
||||
error("solveGaussJordan singular matrix");
|
||||
double pivinv = 1 / A[icol][icol];
|
||||
A[icol][icol] = 1;
|
||||
for (int j = 0; j < n; j++)
|
||||
A[icol][j] *= pivinv;
|
||||
for (int j = 0; j < m; j++)
|
||||
B[icol][j] *= pivinv;
|
||||
|
||||
for (int j = 0; j < n; j++) {
|
||||
if (j != icol) {
|
||||
double mult = A[j][icol];
|
||||
for (int k = 0; k < n; k++)
|
||||
A[j][k] -= A[icol][k] * mult;
|
||||
for (int k = 0; k < m; k++)
|
||||
B[j][k] -= B[icol][k] * mult;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
for (int i = n - 1; i >= 0; i--) {
|
||||
if (indxr[i] != indxc[i]) {
|
||||
for (int j = 0; j < n; j++)
|
||||
std::swap(A[j][indxr[i]], A[j][indxc[i]]);
|
||||
}
|
||||
}
|
||||
|
||||
delete [] indxc;
|
||||
delete [] indxr;
|
||||
delete [] ipiv;
|
||||
}
|
||||
|
||||
void ludcmp(Matrix& A, int *indx, double& d) {
|
||||
const double TINY=1.0e-20;
|
||||
int i,imax,j,k;
|
||||
double big,dum,sum,temp;
|
||||
|
||||
int n=A.numRows();
|
||||
Vector vv(n);
|
||||
d=1.0;
|
||||
for (i=0;i<n;i++) {
|
||||
big=0.0;
|
||||
for (j=0;j<n;j++)
|
||||
if ((temp=std::abs(A[i][j])) > big) big=temp;
|
||||
if (big == 0.0) error("Singular matrix in routine ludcmp");
|
||||
vv[i]=1.0/big;
|
||||
}
|
||||
for (j=0;j<n;j++) {
|
||||
for (i=0;i<j;i++) {
|
||||
sum=A[i][j];
|
||||
for (k=0;k<i;k++) sum -= A[i][k]*A[k][j];
|
||||
A[i][j]=sum;
|
||||
}
|
||||
big=0.0;
|
||||
for (i=j;i<n;i++) {
|
||||
sum=A[i][j];
|
||||
for (k=0;k<j;k++) sum -= A[i][k]*A[k][j];
|
||||
A[i][j]=sum;
|
||||
if ((dum=vv[i]*std::abs(sum)) >= big) {
|
||||
big=dum;
|
||||
imax=i;
|
||||
}
|
||||
}
|
||||
if (j != imax) {
|
||||
for (k=0;k<n;k++) {
|
||||
dum=A[imax][k];
|
||||
A[imax][k]=A[j][k];
|
||||
A[j][k]=dum;
|
||||
}
|
||||
d = -d;
|
||||
vv[imax]=vv[j];
|
||||
}
|
||||
indx[j]=imax;
|
||||
if (A[j][j] == 0.0) A[j][j]=TINY;
|
||||
if (j != n-1) {
|
||||
dum=1.0/(A[j][j]);
|
||||
for (i=j+1;i<n;i++) A[i][j] *= dum;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void lubksb(const Matrix& A, int *indx, Matrix& B) {
|
||||
int i,ii=0,ip,j;
|
||||
double sum;
|
||||
|
||||
int n=A.numRows();
|
||||
for (int k = 0; k < B.numCols(); k++) {
|
||||
for (i=0;i<n;i++) {
|
||||
ip=indx[i];
|
||||
sum=B[ip][k];
|
||||
B[ip][k]=B[i][k];
|
||||
if (ii != 0)
|
||||
for (j=ii-1;j<i;j++) sum -= A[i][j]*B[j][k];
|
||||
else if (sum != 0.0)
|
||||
ii=i+1;
|
||||
B[i][k]=sum;
|
||||
}
|
||||
for (i=n-1;i>=0;i--) {
|
||||
sum=B[i][k];
|
||||
for (j=i+1;j<n;j++) sum -= A[i][j]*B[j][k];
|
||||
B[i][k]=sum/A[i][i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void solveLUDecompose(Matrix& A, Matrix& B) {
|
||||
|
||||
int n = A.numRows();
|
||||
if (A.numCols() != n)
|
||||
error("matrix dimension mismatch");
|
||||
if (B.numRows() != n)
|
||||
error("matrix dimension mismatch");
|
||||
int m = B.numCols();
|
||||
|
||||
int *indx = new int [n];
|
||||
double d;
|
||||
ludcmp(A, indx, d);
|
||||
lubksb(A, indx, B);
|
||||
|
||||
delete [] indx;
|
||||
}
|
||||
|
||||
void reduceHouseholder(Matrix& A, Vector& d, Vector& e) {
|
||||
|
||||
int n = d.dimension();
|
||||
|
||||
for (int i = n - 1; i > 0; i--) {
|
||||
int l = i - 1;
|
||||
double h = 0;
|
||||
double scale = 0;
|
||||
if (l > 0) {
|
||||
for (int k = 0; k <= l; k++)
|
||||
scale += std::abs(A[i][k]);
|
||||
if (scale == 0.0)
|
||||
e[i] = A[i][l];
|
||||
else {
|
||||
for (int k = 0; k <= l; k++) {
|
||||
A[i][k] /= scale;
|
||||
h += A[i][k] * A[i][k];
|
||||
}
|
||||
double f = A[i][l];
|
||||
double g = (f >= 0.0 ? -std::sqrt(h) : std::sqrt(h));
|
||||
e[i] = scale * g;
|
||||
h -= f * g;
|
||||
A[i][l] = f - g;
|
||||
f = 0.0;
|
||||
for (int j = 0; j <= l; j++) {
|
||||
A[j][i] = A[i][j] / h;
|
||||
g = 0.0;
|
||||
for (int k = 0; k <= j; k++)
|
||||
g += A[j][k] * A[i][k];
|
||||
for (int k = j + 1; k <= l; k++)
|
||||
g += A[k][j] * A[i][k];
|
||||
e[j] = g / h;
|
||||
f += e[j] * A[i][j];
|
||||
}
|
||||
double hh = f / (h + h);
|
||||
for (int j = 0; j <= l; j++) {
|
||||
f = A[i][j];
|
||||
e[j] = g = e[j] - hh * f;
|
||||
for (int k = 0; k <= j; k++)
|
||||
A[j][k] -= f * e[k] + g * A[i][k];
|
||||
}
|
||||
}
|
||||
} else
|
||||
e[i] = A[i][l];
|
||||
d[i] = h;
|
||||
}
|
||||
d[0] = 0.0;
|
||||
e[0]=0.0;
|
||||
for (int i = 0; i < n; i++) {
|
||||
if (d[i] != 0.0) {
|
||||
for (int j = 0; j < i; j++) {
|
||||
double g = 0;
|
||||
for (int k = 0; k < i; k++)
|
||||
g += A[i][k] * A[k][j];
|
||||
for (int k = 0; k < i; k++)
|
||||
A[k][j] -= g * A[k][i];
|
||||
}
|
||||
}
|
||||
d[i] = A[i][i];
|
||||
A[i][i] = 1.0;
|
||||
for (int j = 0; j < i; j++)
|
||||
A[j][i] = A[i][j] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
static double pythag(double a, double b) {
|
||||
double absa = std::abs(a);
|
||||
double absb = std::abs(b);
|
||||
if (absa > absb) {
|
||||
double ratio = absb / absa;
|
||||
return absa * std::sqrt(1 + ratio * ratio);
|
||||
} else {
|
||||
if (absb == 0.0)
|
||||
return 0.0;
|
||||
else {
|
||||
double ratio = absa / absb;
|
||||
return absb * std::sqrt(1 + ratio * ratio);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
inline double sign(double a, double b) {
|
||||
if (b >= 0.0) {
|
||||
if (a >= 0.0)
|
||||
return a;
|
||||
else
|
||||
return -a;
|
||||
} else {
|
||||
if (a >= 0.0)
|
||||
return -a;
|
||||
else
|
||||
return a;
|
||||
}
|
||||
}
|
||||
|
||||
void solveTQLI(Vector& d, Vector& e, Matrix& z) {
|
||||
|
||||
int n = d.dimension();
|
||||
for (int i = 1; i < n; i++)
|
||||
e[i-1] = e[i];
|
||||
e[n-1] = 0.0;
|
||||
for (int l = 0; l < n; l++) {
|
||||
int iter = 0, m;
|
||||
do {
|
||||
for (m = l ; m < n-1; m++) {
|
||||
double dd = std::abs(d[m]) + std::abs(d[m+1]);
|
||||
if ((std::abs(e[m]) + dd) == dd)
|
||||
break;
|
||||
}
|
||||
if (m != l) {
|
||||
if (iter++ == 30)
|
||||
error("Too many iterations in solveTQLI");
|
||||
double g = (d[l+1] - d[l]) / (2.0 * e[l]);
|
||||
double r = pythag(g, 1.0);
|
||||
g = d[m] - d[l] + e[l] / (g + sign(r, g));
|
||||
double s = 1.0;
|
||||
double c = 1.0;
|
||||
double p = 0.0;
|
||||
int i;
|
||||
for (i = m-1; i >= l; i--) {
|
||||
double f = s * e[i];
|
||||
double b = c * e[i];
|
||||
e[i+1] = r = pythag(f, g);
|
||||
if (r == 0.0) {
|
||||
d[i+1] -= p;
|
||||
e[m] = 0.0;
|
||||
break;
|
||||
}
|
||||
s = f / r;
|
||||
c = g / r;
|
||||
g = d[i+1] - p;
|
||||
r = (d[i] - g) * s + 2.0 * c * b;
|
||||
d[i+1] = g + (p = s * r);
|
||||
g = c * r - b;
|
||||
for (int k = 0; k < n; k++) {
|
||||
f = z[k][i+1];
|
||||
z[k][i+1] = s * z[k][i] + c * f;
|
||||
z[k][i] = c * z[k][i] - s * f;
|
||||
}
|
||||
}
|
||||
if (r == 0.0 && i >= l)
|
||||
continue;
|
||||
d[l] -= p;
|
||||
e[l] = g;
|
||||
e[m] = 0.0;
|
||||
}
|
||||
} while (m != l);
|
||||
}
|
||||
}
|
||||
|
||||
void sortEigen(Vector& d, Matrix& z) {
|
||||
|
||||
// sorts eigenvalues and eigenvector in descending order
|
||||
int n = d.dimension();
|
||||
if (z.numRows() != n || z.numCols() != n)
|
||||
error("Bad vector, matrix dimensions in sortEigen");
|
||||
for (int i = 0; i < n - 1; i++) {
|
||||
int k = i;
|
||||
double p = d[k];
|
||||
for (int j = i; j < n; j++)
|
||||
if (d[j] >= p)
|
||||
p = d[k = j];
|
||||
if (k != i) {
|
||||
d[k] = d[i];
|
||||
d[i] = p;
|
||||
for (int j = 0; j < n; j++) {
|
||||
p = z[j][i];
|
||||
z[j][i] = z[j][k];
|
||||
z[j][k] = p;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void singularValueDecompose(Matrix& a, Vector& w, Matrix& v) {
|
||||
|
||||
bool flag;
|
||||
int i,its,j,jj,k,l,nm;
|
||||
double anorm,c,f,g,h,s,scale,x,y,z;
|
||||
|
||||
int m=a.numRows();
|
||||
int n=a.numCols();
|
||||
Vector rv1(n);
|
||||
g=scale=anorm=0.0;
|
||||
for (i=0;i<n;i++) {
|
||||
l=i+2;
|
||||
rv1[i]=scale*g;
|
||||
g=s=scale=0.0;
|
||||
if (i < m) {
|
||||
for (k=i;k<m;k++) scale += std::abs(a[k][i]);
|
||||
if (scale != 0.0) {
|
||||
for (k=i;k<m;k++) {
|
||||
a[k][i] /= scale;
|
||||
s += a[k][i]*a[k][i];
|
||||
}
|
||||
f=a[i][i];
|
||||
g = -sign(std::sqrt(s),f);
|
||||
h=f*g-s;
|
||||
a[i][i]=f-g;
|
||||
for (j=l-1;j<n;j++) {
|
||||
for (s=0.0,k=i;k<m;k++) s += a[k][i]*a[k][j];
|
||||
f=s/h;
|
||||
for (k=i;k<m;k++) a[k][j] += f*a[k][i];
|
||||
}
|
||||
for (k=i;k<m;k++) a[k][i] *= scale;
|
||||
}
|
||||
}
|
||||
w[i]=scale *g;
|
||||
g=s=scale=0.0;
|
||||
if (i+1 <= m && i != n) {
|
||||
for (k=l-1;k<n;k++) scale += std::abs(a[i][k]);
|
||||
if (scale != 0.0) {
|
||||
for (k=l-1;k<n;k++) {
|
||||
a[i][k] /= scale;
|
||||
s += a[i][k]*a[i][k];
|
||||
}
|
||||
f=a[i][l-1];
|
||||
g = -sign(std::sqrt(s),f);
|
||||
h=f*g-s;
|
||||
a[i][l-1]=f-g;
|
||||
for (k=l-1;k<n;k++) rv1[k]=a[i][k]/h;
|
||||
for (j=l-1;j<m;j++) {
|
||||
for (s=0.0,k=l-1;k<n;k++) s += a[j][k]*a[i][k];
|
||||
for (k=l-1;k<n;k++) a[j][k] += s*rv1[k];
|
||||
}
|
||||
for (k=l-1;k<n;k++) a[i][k] *= scale;
|
||||
}
|
||||
}
|
||||
anorm=std::max(anorm,(std::abs(w[i])+std::abs(rv1[i])));
|
||||
}
|
||||
for (i=n-1;i>=0;i--) {
|
||||
if (i < n-1) {
|
||||
if (g != 0.0) {
|
||||
for (j=l;j<n;j++)
|
||||
v[j][i]=(a[i][j]/a[i][l])/g;
|
||||
for (j=l;j<n;j++) {
|
||||
for (s=0.0,k=l;k<n;k++) s += a[i][k]*v[k][j];
|
||||
for (k=l;k<n;k++) v[k][j] += s*v[k][i];
|
||||
}
|
||||
}
|
||||
for (j=l;j<n;j++) v[i][j]=v[j][i]=0.0;
|
||||
}
|
||||
v[i][i]=1.0;
|
||||
g=rv1[i];
|
||||
l=i;
|
||||
}
|
||||
for (i=std::min(m,n)-1;i>=0;i--) {
|
||||
l=i+1;
|
||||
g=w[i];
|
||||
for (j=l;j<n;j++) a[i][j]=0.0;
|
||||
if (g != 0.0) {
|
||||
g=1.0/g;
|
||||
for (j=l;j<n;j++) {
|
||||
for (s=0.0,k=l;k<m;k++) s += a[k][i]*a[k][j];
|
||||
f=(s/a[i][i])*g;
|
||||
for (k=i;k<m;k++) a[k][j] += f*a[k][i];
|
||||
}
|
||||
for (j=i;j<m;j++) a[j][i] *= g;
|
||||
} else for (j=i;j<m;j++) a[j][i]=0.0;
|
||||
++a[i][i];
|
||||
}
|
||||
for (k=n-1;k>=0;k--) {
|
||||
for (its=0;its<30;its++) {
|
||||
flag=true;
|
||||
for (l=k;l>=0;l--) {
|
||||
nm=l-1;
|
||||
if ((std::abs(rv1[l])+anorm) == anorm) {
|
||||
flag=false;
|
||||
break;
|
||||
}
|
||||
if ((std::abs(w[nm])+anorm) == anorm) break;
|
||||
}
|
||||
if (flag) {
|
||||
c=0.0;
|
||||
s=1.0;
|
||||
for (i=l-1;i<k+1;i++) {
|
||||
f=s*rv1[i];
|
||||
rv1[i]=c*rv1[i];
|
||||
if ((std::abs(f)+anorm) == anorm) break;
|
||||
g=w[i];
|
||||
h=pythag(f,g);
|
||||
w[i]=h;
|
||||
h=1.0/h;
|
||||
c=g*h;
|
||||
s = -f*h;
|
||||
for (j=0;j<m;j++) {
|
||||
y=a[j][nm];
|
||||
z=a[j][i];
|
||||
a[j][nm]=y*c+z*s;
|
||||
a[j][i]=z*c-y*s;
|
||||
}
|
||||
}
|
||||
}
|
||||
z=w[k];
|
||||
if (l == k) {
|
||||
if (z < 0.0) {
|
||||
w[k] = -z;
|
||||
for (j=0;j<n;j++) v[j][k] = -v[j][k];
|
||||
}
|
||||
break;
|
||||
}
|
||||
if (its == 30) error("singularValueDecomposition: no convergence in 30 iterations");
|
||||
x=w[l];
|
||||
nm=k-1;
|
||||
y=w[nm];
|
||||
g=rv1[nm];
|
||||
h=rv1[k];
|
||||
f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
|
||||
g=pythag(f,1.0);
|
||||
f=((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x;
|
||||
c=s=1.0;
|
||||
for (j=l;j<=nm;j++) {
|
||||
i=j+1;
|
||||
g=rv1[i];
|
||||
y=w[i];
|
||||
h=s*g;
|
||||
g=c*g;
|
||||
z=pythag(f,h);
|
||||
rv1[j]=z;
|
||||
c=f/z;
|
||||
s=h/z;
|
||||
f=x*c+g*s;
|
||||
g = g*c-x*s;
|
||||
h=y*s;
|
||||
y *= c;
|
||||
for (jj=0;jj<n;jj++) {
|
||||
x=v[jj][j];
|
||||
z=v[jj][i];
|
||||
v[jj][j]=x*c+z*s;
|
||||
v[jj][i]=z*c-x*s;
|
||||
}
|
||||
z=pythag(f,h);
|
||||
w[j]=z;
|
||||
if (z != 0.0) {
|
||||
z=1.0/z;
|
||||
c=f*z;
|
||||
s=h*z;
|
||||
}
|
||||
f=c*g+s*y;
|
||||
x=c*y-s*g;
|
||||
for (jj=0;jj<m;jj++) {
|
||||
y=a[jj][j];
|
||||
z=a[jj][i];
|
||||
a[jj][j]=y*c+z*s;
|
||||
a[jj][i]=z*c-y*s;
|
||||
}
|
||||
}
|
||||
rv1[l]=0.0;
|
||||
rv1[k]=f;
|
||||
w[k]=x;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Vector solveEigenSymmetric(Matrix& A) {
|
||||
int n = A.numRows();
|
||||
Vector d(n), e(n);
|
||||
reduceHouseholder(A, d, e);
|
||||
solveTQLI(d, e, A);
|
||||
sortEigen(d, A);
|
||||
return d;
|
||||
}
|
||||
|
||||
void lineSearch(Vector& xold, double fold, Vector& g, Vector& p, Vector& x,
|
||||
double& f, double stpmax, bool& check, double (*func)(Vector&))
|
||||
{
|
||||
const double ALF=1.0e-4, TOLX=std::numeric_limits<double>::epsilon();
|
||||
int i;
|
||||
double a,alam,alam2=0.0,alamin,b,disc,f2=0.0;
|
||||
double rhs1,rhs2,slope,sum,temp,test,tmplam;
|
||||
|
||||
int n = xold.dimension();
|
||||
check = false;
|
||||
sum=0.0;
|
||||
for (i=0;i<n;i++) sum += p[i]*p[i];
|
||||
sum=std::sqrt(sum);
|
||||
if (sum > stpmax)
|
||||
for (i=0;i<n;i++) p[i] *= stpmax/sum;
|
||||
slope = 0.0;
|
||||
for (i=0;i<n;i++)
|
||||
slope += g[i]*p[i];
|
||||
if (slope >= 0.0) error("Roundoff problem in lineSearch");
|
||||
test=0.0;
|
||||
for (i=0;i<n;i++) {
|
||||
temp=std::abs(p[i])/std::max(std::abs(xold[i]),1.0);
|
||||
if (temp > test) test=temp;
|
||||
}
|
||||
alamin=TOLX/test;
|
||||
alam=1.0;
|
||||
for (;;) {
|
||||
for (i=0;i<n;i++) x[i]=xold[i]+alam*p[i];
|
||||
f=func(x);
|
||||
if (alam < alamin) {
|
||||
for (i=0;i<n;i++) x[i]=xold[i];
|
||||
check = true;
|
||||
return;
|
||||
} else if (f <= fold+ALF*alam*slope) return;
|
||||
else {
|
||||
if (alam == 1.0)
|
||||
tmplam = -slope/(2.0*(f-fold-slope));
|
||||
else {
|
||||
rhs1 = f-fold-alam*slope;
|
||||
rhs2=f2-fold-alam2*slope;
|
||||
a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
|
||||
b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
|
||||
if (a == 0.0) tmplam = -slope/(2.0*b);
|
||||
else {
|
||||
disc=b*b-3.0*a*slope;
|
||||
if (disc < 0.0) tmplam=0.5*alam;
|
||||
else if (b <= 0.0) tmplam=(-b+std::sqrt(disc))/(3.0*a);
|
||||
else tmplam=-slope/(b+std::sqrt(disc));
|
||||
}
|
||||
if (tmplam > 0.5*alam)
|
||||
tmplam=0.5*alam;
|
||||
}
|
||||
}
|
||||
alam2=alam;
|
||||
f2 = f;
|
||||
alam=std::max(tmplam,0.1*alam);
|
||||
}
|
||||
}
|
||||
|
||||
void minimizeBFGS(Vector& p, const double gtol, int& iter, double& fret,
|
||||
double (*func)(Vector&), void (*dfunc)(Vector&, Vector&)) {
|
||||
const int ITMAX = 200;
|
||||
const double EPS = std::numeric_limits<double>::epsilon();
|
||||
const double TOLX = 4*EPS, STPMX = 100.0;
|
||||
bool check;
|
||||
int i,its,j;
|
||||
double den,fac,fad,fae,fp,stpmax,sum=0.0,sumdg,sumxi,temp,test;
|
||||
|
||||
int n = p.dimension();
|
||||
Vector dg(n),g(n),hdg(n),pnew(n),xi(n);
|
||||
Matrix hessin(n, n);
|
||||
|
||||
fp=func(p);
|
||||
dfunc(p,g);
|
||||
for (i=0;i<n;i++) {
|
||||
for (j=0;j<n;j++) hessin[i][j]=0.0;
|
||||
hessin[i][i]=1.0;
|
||||
xi[i] = -g[i];
|
||||
sum += p[i]*p[i];
|
||||
}
|
||||
stpmax=STPMX*std::max(std::sqrt(sum),double(n));
|
||||
for (its=0;its<ITMAX;its++) {
|
||||
iter=its;
|
||||
lineSearch(p,fp,g,xi,pnew,fret,stpmax,check,func);
|
||||
fp = fret;
|
||||
for (i=0;i<n;i++) {
|
||||
xi[i]=pnew[i]-p[i];
|
||||
p[i]=pnew[i];
|
||||
}
|
||||
test=0.0;
|
||||
for (i=0;i<n;i++) {
|
||||
temp=std::abs(xi[i])/std::max(std::abs(p[i]),1.0);
|
||||
if (temp > test) test=temp;
|
||||
}
|
||||
if (test < TOLX)
|
||||
return;
|
||||
for (i=0;i<n;i++) dg[i]=g[i];
|
||||
dfunc(p,g);
|
||||
test=0.0;
|
||||
den=std::max(fret,1.0);
|
||||
for (i=0;i<n;i++) {
|
||||
temp=std::abs(g[i])*std::max(std::abs(p[i]),1.0)/den;
|
||||
if (temp > test) test=temp;
|
||||
}
|
||||
if (test < gtol)
|
||||
return;
|
||||
for (i=0;i<n;i++) dg[i]=g[i]-dg[i];
|
||||
for (i=0;i<n;i++) {
|
||||
hdg[i]=0.0;
|
||||
for (j=0;j<n;j++) hdg[i] += hessin[i][j]*dg[j];
|
||||
}
|
||||
fac=fae=sumdg=sumxi=0.0;
|
||||
for (i=0;i<n;i++) {
|
||||
fac += dg[i]*xi[i];
|
||||
fae += dg[i]*hdg[i];
|
||||
sumdg += dg[i]*dg[i];
|
||||
sumxi += xi[i]*xi[i];
|
||||
}
|
||||
if (fac > std::sqrt(EPS*sumdg*sumxi)) {
|
||||
fac=1.0/fac;
|
||||
fad=1.0/fae;
|
||||
for (i=0;i<n;i++) dg[i]=fac*xi[i]-fad*hdg[i];
|
||||
for (i=0;i<n;i++) {
|
||||
for (j=i;j<n;j++) {
|
||||
hessin[i][j] += fac*xi[i]*xi[j]
|
||||
-fad*hdg[i]*hdg[j]+fae*dg[i]*dg[j];
|
||||
hessin[j][i]=hessin[i][j];
|
||||
}
|
||||
}
|
||||
}
|
||||
for (i=0;i<n;i++) {
|
||||
xi[i]=0.0;
|
||||
for (j=0;j<n;j++) xi[i] -= hessin[i][j]*g[j];
|
||||
}
|
||||
}
|
||||
error("Too many iterations in minimizeBFGS");
|
||||
}
|
||||
|
||||
} /* end namespace cpl */
|
||||
|
|
@ -0,0 +1,86 @@
|
|||
#ifndef CPL_MATRIX_HPP
|
||||
#define CPL_MATRIX_HPP
|
||||
|
||||
namespace cpl {
|
||||
|
||||
class Vector;
|
||||
|
||||
class Matrix {
|
||||
public:
|
||||
|
||||
Matrix(int rows=1, int cols=1, double d=0);
|
||||
|
||||
Matrix(const Matrix& m);
|
||||
|
||||
~Matrix();
|
||||
|
||||
int numRows() const { return rows; }
|
||||
|
||||
int numCols() const { return cols; }
|
||||
|
||||
const double* operator[](const int row) const;
|
||||
|
||||
double* operator[](const int row);
|
||||
|
||||
double& operator()(int row, int col);
|
||||
|
||||
Matrix& operator=(const Matrix& m);
|
||||
|
||||
Matrix& operator=(const double d);
|
||||
|
||||
Matrix& operator += (const Matrix& m);
|
||||
|
||||
Matrix& operator -= (const Matrix& m);
|
||||
|
||||
Matrix& operator *= (double d);
|
||||
|
||||
Matrix& operator /= (double d);
|
||||
|
||||
Matrix transpose();
|
||||
|
||||
friend std::ostream& operator<<(std::ostream& os, const Matrix& mat);
|
||||
|
||||
private:
|
||||
int rows;
|
||||
int cols;
|
||||
double **m;
|
||||
};
|
||||
|
||||
inline Matrix operator + (const Matrix& m) {
|
||||
return m;
|
||||
}
|
||||
|
||||
extern Matrix operator - (const Matrix& m);
|
||||
|
||||
extern Matrix operator * (const Matrix& m, double d);
|
||||
|
||||
extern Matrix operator * (double d, const Matrix& m);
|
||||
|
||||
extern Matrix operator / (const Matrix& m, double d);
|
||||
|
||||
extern Matrix operator + (const Matrix& m1, const Matrix& m2);
|
||||
|
||||
extern Matrix operator - (const Matrix& m1, const Matrix& m2);
|
||||
|
||||
extern Matrix operator * (const Matrix& m1, const Matrix& m2);
|
||||
|
||||
extern void solveGaussJordan(Matrix& A, Matrix& B);
|
||||
|
||||
extern void solveLUDecompose(Matrix& A, Matrix& B);
|
||||
|
||||
extern void reduceHouseholder(Matrix& A, Vector& d, Vector& e);
|
||||
|
||||
extern void solveTQLI(Vector& d, Vector&e, Matrix& z);
|
||||
|
||||
extern void sortEigen(Vector& d, Matrix& z);
|
||||
|
||||
extern Vector solveEigenSymmetric(Matrix& A);
|
||||
|
||||
extern void singularValueDecompose(Matrix& a, Vector& w, Matrix& v);
|
||||
|
||||
extern void minimizeBFGS(Vector& p, const double gtol, int& iter, double& fret,
|
||||
double (*func)(Vector&), void (*dfunc)(Vector&, Vector&));
|
||||
|
||||
} /* end namespace cpl */
|
||||
|
||||
#endif /* CPL_MATRIX_HPP */
|
||||
|
|
@ -0,0 +1,136 @@
|
|||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <cstdlib>
|
||||
#include <iostream>
|
||||
using namespace std;
|
||||
|
||||
#include "vector.hpp"
|
||||
#include "odeint.hpp"
|
||||
|
||||
namespace cpl {
|
||||
|
||||
// Fourth order Runge-Kutta
|
||||
void RK4Step(Vector& x, double tau,
|
||||
Vector derivs(const Vector&))
|
||||
{
|
||||
Vector k1 = tau * derivs(x);
|
||||
Vector k2 = tau * derivs(x + 0.5 * k1);
|
||||
Vector k3 = tau * derivs(x + 0.5 * k2);
|
||||
Vector k4 = tau * derivs(x + k3);
|
||||
x += (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;
|
||||
}
|
||||
|
||||
// Adaptive step size control using Runge-Kutta and step doubling
|
||||
void adaptiveRK4Step(Vector& x, double& tau, double accuracy,
|
||||
Vector derivs(const Vector&))
|
||||
{
|
||||
const double SAFETY = 0.9, PGROW = -0.2, PSHRINK = -0.25,
|
||||
ERRCON = 1.89E-4, TINY = 1.0e-30;
|
||||
int n = x.dimension();
|
||||
Vector x_half(n), x_full(n), Delta(n);
|
||||
Vector scale = derivs(x);
|
||||
for (int i = 0; i < n; i++)
|
||||
scale[i] = abs(x[i]) + abs(scale[i] * tau) + TINY;
|
||||
double err_max;
|
||||
while (true) {
|
||||
// take two half steps
|
||||
double tau_half = tau / 2;
|
||||
x_half = x;
|
||||
RK4Step(x_half, tau_half, derivs);
|
||||
RK4Step(x_half, tau_half, derivs);
|
||||
// take full step
|
||||
x_full = x;
|
||||
RK4Step(x_full, tau, derivs);
|
||||
// estimate error
|
||||
Delta = x_half - x_full;
|
||||
err_max = 0;
|
||||
for (int i = 0; i < n; i++)
|
||||
err_max = max(err_max, abs(Delta[i]) / scale[i]);
|
||||
err_max /= accuracy;
|
||||
if (err_max <= 1.0)
|
||||
break;
|
||||
double tau_temp = SAFETY * tau * pow(err_max, PSHRINK);
|
||||
if (tau >= 0.0)
|
||||
tau = max(tau_temp, 0.1 * tau);
|
||||
else
|
||||
tau = min(tau_temp, 0.1 * tau);
|
||||
if (abs(tau) == 0.0) {
|
||||
cerr << "adaptiveRK4Step: step size underflow\naborting ..."
|
||||
<< endl;
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
}
|
||||
tau *= (err_max > ERRCON ? SAFETY * pow(err_max, PGROW) : 5.0);
|
||||
x = x_half + Delta / 15.0;
|
||||
}
|
||||
|
||||
// Runge-Kutta-Cash-Karp including error estimate
|
||||
static void rkck(Vector& x, double tau,
|
||||
Vector derivs(const Vector&), Vector& x_err)
|
||||
{
|
||||
const double b21 = 1.0/5.0, b31 = 3.0/40.0, b41 = 3.0/10.0,
|
||||
b51 = -11.0/54.0, b61 = 1631.0/55296.0, b32 = 9.0/40.0,
|
||||
b42 = -9.0/10.0, b52 = 5.0/2.0, b62 = 175.0/512.0, b43 = 6.0/5.0,
|
||||
b53 = -70.0/27.0, b63 = 575.0/13824.0, b54 = 35.0/27.0,
|
||||
b64 = 44275.0/110592.0, b65 = 253.0/4096.0, c1 = 37.0/378.0,
|
||||
c2 = 0.0, c3 = 250.0/621.0, c4 = 125.0/594.0, c5 = 0.0,
|
||||
c6 = 512.0/1771.0, dc1 = c1 - 2825.0/27648.0, dc2 = c2 - 0.0,
|
||||
dc3 = c3 - 18575.0/48384.0, dc4 = c4 - 13525.0/55296.0,
|
||||
dc5 = c5 - 277.0/14336.0, dc6 = c6 - 1.0/4.0;
|
||||
|
||||
Vector
|
||||
k1 = tau * derivs(x),
|
||||
k2 = tau * derivs(x + b21*k1),
|
||||
k3 = tau * derivs(x + b31*k1 + b32*k2),
|
||||
k4 = tau * derivs(x + b41*k1 + b42*k2 + b43*k3),
|
||||
k5 = tau * derivs(x + b51*k1 + b52*k2 + b53*k3 + b54*k4),
|
||||
k6 = tau * derivs(x + b61*k1 + b62*k2 + b63*k3 + b64*k4 + b65*k5);
|
||||
x += c1*k1 + c2*k2 + c3*k3 + c4*k4 + c5*k5 + c6*k6;
|
||||
x_err = dc1*k1 + dc2*k2 + dc3*k3 + dc4*k4 + dc5*k5 + dc6*k6;
|
||||
}
|
||||
|
||||
// Runge-Kutta-Cash-Karp step
|
||||
void RKCKStep(Vector& x, double tau, Vector derivs(const Vector&))
|
||||
{
|
||||
Vector x_err(x.dimension());
|
||||
rkck(x, tau, derivs, x_err);
|
||||
}
|
||||
|
||||
// Adaptive step size control using Runge-Kutta-Cash-Karp
|
||||
void adaptiveRKCKStep(Vector& x, double& tau, double accuracy,
|
||||
Vector derivs(const Vector&))
|
||||
{
|
||||
const double SAFETY = 0.9, PGROW = -0.2, PSHRINK = -0.25,
|
||||
ERRCON = 1.89E-4, TINY = 1.0e-30;
|
||||
int n = x.dimension();
|
||||
Vector x_err(n), x_temp(n);
|
||||
Vector scale = derivs(x);
|
||||
for (int i = 0; i < n; i++)
|
||||
scale[i] = abs(x[i]) + abs(scale[i] * tau) + TINY;
|
||||
double err_max;
|
||||
while (true) {
|
||||
// take Cash-Karp step including error estimate
|
||||
x_temp = x;
|
||||
rkck(x_temp, tau, derivs, x_err);
|
||||
err_max = 0;
|
||||
for (int i = 0; i < n; i++)
|
||||
err_max = max(err_max, abs(x_err[i]) / scale[i]);
|
||||
err_max /= accuracy;
|
||||
if (err_max <= 1.0)
|
||||
break;
|
||||
double tau_temp = SAFETY * tau * pow(err_max, PSHRINK);
|
||||
if (tau >= 0.0)
|
||||
tau = max(tau_temp, 0.1 * tau);
|
||||
else
|
||||
tau = min(tau_temp, 0.1 * tau);
|
||||
if (abs(tau) == 0.0) {
|
||||
cerr << "adaptiveRKCKStep: step size underflow\naborting ..."
|
||||
<< endl;
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
}
|
||||
tau *= (err_max > ERRCON ? SAFETY * pow(err_max, PGROW) : 5.0);
|
||||
x = x_temp;
|
||||
}
|
||||
|
||||
} /* end namespace cpl */
|
||||
|
|
@ -0,0 +1,44 @@
|
|||
#ifndef CPL_ODEINT_HPP
|
||||
#define CPL_ODEINT_HPP
|
||||
|
||||
#include <iostream>
|
||||
|
||||
namespace cpl {
|
||||
|
||||
class Vector;
|
||||
|
||||
// ODE integration routines adapted from Numerical Recipes
|
||||
|
||||
// take a single 4th order Runge-Kutta step
|
||||
extern void RK4Step(
|
||||
Vector& x, // extended solution vector
|
||||
double tau, // step size
|
||||
Vector derivs(const Vector& x) // extended derivative vector
|
||||
);
|
||||
|
||||
// take one adaptive RK4 step using step doubling
|
||||
extern void adaptiveRK4Step(
|
||||
Vector& x, // extended solution vector
|
||||
double& tau, // recommended step size
|
||||
double accuracy, // desired solution accuracy
|
||||
Vector derivs(const Vector&) // extended derivative vector
|
||||
);
|
||||
|
||||
// take one Runge-Kutta-Cash-Karp step
|
||||
extern void RKCKStep(
|
||||
Vector& x, // extended solution vector
|
||||
double tau, // step size
|
||||
Vector derivs(const Vector& x) // derivative vector
|
||||
);
|
||||
|
||||
// take one adaptive RKCK step using embedded error estimate
|
||||
extern void adaptiveRKCKStep(
|
||||
Vector& x, // extended solution vector
|
||||
double& tau, // recommended step size
|
||||
double accuracy, // desired solution accuracy
|
||||
Vector derivs(const Vector&) // extended derivative vector
|
||||
);
|
||||
|
||||
} /* end namespace cpl */
|
||||
|
||||
#endif /* CPL_ODEINT_HPP */
|
||||
|
|
@ -0,0 +1,134 @@
|
|||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <cstdlib>
|
||||
#include <iostream>
|
||||
using namespace std;
|
||||
|
||||
#include "optimize.hpp"
|
||||
|
||||
static int sign_of_func;
|
||||
static double (*user_func)(double);
|
||||
static double func(double x) {
|
||||
return sign_of_func * (*user_func)(x);
|
||||
}
|
||||
|
||||
inline void shft3(double &a, double &b, double &c, const double d) {
|
||||
a=b; b=c; c=d;
|
||||
}
|
||||
|
||||
inline double SIGN(const double& a, const double& b) {
|
||||
return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);
|
||||
}
|
||||
|
||||
namespace cpl {
|
||||
|
||||
double Optimizer::findMaximum(const double a, const double b,
|
||||
double f(double), double tol, double& xmin) {
|
||||
sign_of_func = -1;
|
||||
user_func = f;
|
||||
double ax = a, bx = b, cx, fa, fb, fc;
|
||||
mnbrak(ax, bx, cx, fa, fb, fc, func);
|
||||
return -golden(ax, bx, cx, func, tol, xmin);
|
||||
}
|
||||
|
||||
double Optimizer::findMinimum(const double a, const double b,
|
||||
double f(double), double tol, double& xmin) {
|
||||
double ax = a, bx = b, cx, fa, fb, fc;
|
||||
mnbrak(ax, bx, cx, fa, fb, fc, f);
|
||||
return golden(ax, bx, cx, f, tol, xmin);
|
||||
}
|
||||
|
||||
void Optimizer::mnbrak(double &ax, double &bx, double &cx,
|
||||
double &fa, double &fb, double &fc,
|
||||
double func(const double)) {
|
||||
const double GOLD=1.618034,GLIMIT=100.0,TINY=1.0e-20;
|
||||
double ulim,u,r,q,fu;
|
||||
|
||||
fa=func(ax);
|
||||
fb=func(bx);
|
||||
if (fb > fa) {
|
||||
swap(ax,bx);
|
||||
swap(fb,fa);
|
||||
}
|
||||
cx=bx+GOLD*(bx-ax);
|
||||
fc=func(cx);
|
||||
while (fb > fc) {
|
||||
r=(bx-ax)*(fb-fc);
|
||||
q=(bx-cx)*(fb-fa);
|
||||
u=bx-((bx-cx)*q-(bx-ax)*r)/
|
||||
(2.0*SIGN(max(abs(q-r),TINY),q-r));
|
||||
ulim=bx+GLIMIT*(cx-bx);
|
||||
if ((bx-u)*(u-cx) > 0.0) {
|
||||
fu=func(u);
|
||||
if (fu < fc) {
|
||||
ax=bx;
|
||||
bx=u;
|
||||
fa=fb;
|
||||
fb=fu;
|
||||
return;
|
||||
} else if (fu > fb) {
|
||||
cx=u;
|
||||
fc=fu;
|
||||
return;
|
||||
}
|
||||
u=cx+GOLD*(cx-bx);
|
||||
fu=func(u);
|
||||
} else if ((cx-u)*(u-ulim) > 0.0) {
|
||||
fu=func(u);
|
||||
if (fu < fc) {
|
||||
shft3(bx,cx,u,cx+GOLD*(cx-bx));
|
||||
shft3(fb,fc,fu,func(u));
|
||||
}
|
||||
} else if ((u-ulim)*(ulim-cx) >= 0.0) {
|
||||
u=ulim;
|
||||
fu=func(u);
|
||||
} else {
|
||||
u=cx+GOLD*(cx-bx);
|
||||
fu=func(u);
|
||||
}
|
||||
shft3(ax,bx,cx,u);
|
||||
shft3(fa,fb,fc,fu);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
inline void shft2(double &a, double &b, const double c) {
|
||||
a=b; b=c;
|
||||
}
|
||||
|
||||
double Optimizer::golden(const double ax, const double bx, const double cx,
|
||||
double f(const double), const double tol,
|
||||
double &xmin) {
|
||||
const double R=0.61803399,C=1.0-R;
|
||||
double f1,f2,x0,x1,x2,x3;
|
||||
|
||||
x0=ax;
|
||||
x3=cx;
|
||||
if (abs(cx-bx) > abs(bx-ax)) {
|
||||
x1=bx;
|
||||
x2=bx+C*(cx-bx);
|
||||
} else {
|
||||
x2=bx;
|
||||
x1=bx-C*(bx-ax);
|
||||
}
|
||||
f1=f(x1);
|
||||
f2=f(x2);
|
||||
while (abs(x3-x0) > tol*(abs(x1)+abs(x2))) {
|
||||
if (f2 < f1) {
|
||||
shft3(x0,x1,x2,R*x2+C*x3);
|
||||
shft2(f1,f2,f(x2));
|
||||
} else {
|
||||
shft3(x3,x2,x1,R*x1+C*x0);
|
||||
shft2(f2,f1,f(x1));
|
||||
}
|
||||
}
|
||||
if (f1 < f2) {
|
||||
xmin=x1;
|
||||
return f1;
|
||||
} else {
|
||||
xmin=x2;
|
||||
return f2;
|
||||
}
|
||||
}
|
||||
|
||||
} /* end namespace cpl */
|
||||
|
|
@ -0,0 +1,47 @@
|
|||
#ifndef CPL_OPTIMIZE_HPP
|
||||
#define CPL_OPTIMIZE_HPP
|
||||
|
||||
namespace cpl {
|
||||
|
||||
class Optimizer {
|
||||
|
||||
public:
|
||||
|
||||
double findMaximum( // returns f(x) at local maximum
|
||||
const double xa, // guess for one point near desried maximum
|
||||
const double xb, // guess for second point near maximum
|
||||
double f(double), // name of function to be maximized
|
||||
double accuracy, // desired accuracy in x
|
||||
double& xMax); // value of x at local maximum
|
||||
|
||||
double findMinimum( // returns f(x) at local minimum
|
||||
const double xa, // guess for one point near desried minimum
|
||||
const double xb, // guess for second point near minimum
|
||||
double f(double), // name of function to be minimized
|
||||
double accuracy, // desired accuracy in x
|
||||
double& xMin); // value of x at local minimum
|
||||
|
||||
double golden( // numerical recipes routine
|
||||
const double,
|
||||
const double,
|
||||
const double,
|
||||
double f(double),
|
||||
const double,
|
||||
double&
|
||||
);
|
||||
|
||||
void mnbrak( // numerical recipes bracketing routine
|
||||
double&,
|
||||
double&,
|
||||
double&,
|
||||
double&,
|
||||
double&,
|
||||
double&,
|
||||
double f(double)
|
||||
);
|
||||
|
||||
};
|
||||
|
||||
} /* end namespace cpl */
|
||||
|
||||
#endif /* CPL_OPTIMIZE_HPP */
|
||||
|
|
@ -0,0 +1,73 @@
|
|||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <cstdlib>
|
||||
#include <ctime>
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
using namespace std;
|
||||
|
||||
#include "random.hpp"
|
||||
|
||||
namespace cpl {
|
||||
|
||||
int timeSeed() {
|
||||
time_t t; // usually an unsigned long int
|
||||
time(&t); // get seconds since Jan 1, 1900
|
||||
tm* pt = gmtime(&t); // Universal (Greenwich Mean) Time
|
||||
return pt->tm_year + 70 * (pt->tm_mon + 12 * (pt->tm_hour +
|
||||
23 * (pt->tm_min + 59 * pt->tm_sec)));
|
||||
}
|
||||
|
||||
double rand(int& seed, bool set) {
|
||||
if (set)
|
||||
srand(seed);
|
||||
else
|
||||
seed = std::rand();
|
||||
return seed / (RAND_MAX + 1.0);
|
||||
}
|
||||
|
||||
double ranf(int& seed, bool set) {
|
||||
const int IA = 16807, IC = 2147483647, IQ = 127773, IR = 2836;
|
||||
if (!set) {
|
||||
int h = seed / IQ;
|
||||
int l = seed % IQ;
|
||||
seed = IA * l - IR * h;
|
||||
}
|
||||
if (seed <= 0)
|
||||
seed += IC;
|
||||
return seed / double(IC);
|
||||
}
|
||||
|
||||
double qand(int& seed, bool set) {
|
||||
static unsigned long idum;
|
||||
const double TWO_POWER_32 = 4294967296.0;
|
||||
if (set) {
|
||||
idum = (unsigned long) seed;
|
||||
return idum / TWO_POWER_32;
|
||||
}
|
||||
idum = 1664525L * idum + 1013904223L;
|
||||
seed = int(idum);
|
||||
return idum / TWO_POWER_32;
|
||||
}
|
||||
|
||||
double gasdev(int& seed) {
|
||||
static int iset = 0;
|
||||
static double gset;
|
||||
double fac, rsq, v1, v2;
|
||||
if (iset == 0) {
|
||||
do {
|
||||
v1 = 2.0 * ranf(seed) - 1.0;
|
||||
v2 = 2.0 * ranf(seed) - 1.0;
|
||||
rsq = v1 * v1 + v2 * v2;
|
||||
} while (rsq >= 1.0 || rsq == 0.0);
|
||||
fac = sqrt(-2.0 * log(rsq)/rsq);
|
||||
gset = v1 * fac;
|
||||
iset = 1;
|
||||
return v2 * fac;
|
||||
} else {
|
||||
iset = 0;
|
||||
return gset;
|
||||
}
|
||||
}
|
||||
|
||||
} /* end namespace cpl */
|
||||
|
|
@ -0,0 +1,23 @@
|
|||
#ifndef CPL_RANDOM_HPP
|
||||
#define CPL_RANDOM_HPP
|
||||
|
||||
namespace cpl {
|
||||
|
||||
// Random integer seed based on the current time of day
|
||||
extern int timeSeed();
|
||||
|
||||
// Standard C++ Library generator
|
||||
extern double rand(int& seed, bool set=false);
|
||||
|
||||
// Park-Miller generator with Schrage's algorithm to prevent overflow
|
||||
extern double ranf(int& seed, bool set=false);
|
||||
|
||||
// Quick and dirty generator according to Numerical Recipes
|
||||
extern double qand(int& seed, bool set=false);
|
||||
|
||||
// Gaussian distributed random number based on Box-Muller algorithm
|
||||
extern double gasdev(int& seed);
|
||||
|
||||
} /* end namespace cpl */
|
||||
|
||||
#endif /* CPL_RANDOM_HPP */
|
||||
|
|
@ -0,0 +1,167 @@
|
|||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <cstdarg>
|
||||
#include <cstdlib>
|
||||
#include <iostream>
|
||||
using namespace std;
|
||||
|
||||
#include "vector.hpp"
|
||||
|
||||
namespace cpl {
|
||||
|
||||
Vector::Vector(int dim) {
|
||||
v = new double [this->dim = dim];
|
||||
for (int i = 0; i < dim; i++) v[i] = 0;
|
||||
}
|
||||
|
||||
Vector::Vector(double v0, double v1) {
|
||||
v = new double [dim = 2];
|
||||
v[0] = v0;
|
||||
v[1] = v1;
|
||||
}
|
||||
|
||||
Vector::Vector(double v0, double v1, double v2) {
|
||||
v = new double [dim = 3];
|
||||
v[0] = v0;
|
||||
v[1] = v1;
|
||||
v[2] = v2;
|
||||
}
|
||||
|
||||
Vector::Vector(const Vector& dv) {
|
||||
v = new double [dim = dv.dim];
|
||||
for (int i = 0; i < dim; i++) v[i] = dv.v[i];
|
||||
}
|
||||
|
||||
void Vector::resize(const int dimension) {
|
||||
double *old = new double[dim];
|
||||
for (int i = 0; i < dim; i++) old[i] = v[i];
|
||||
delete [] v;
|
||||
v = new double [dimension];
|
||||
for (int i = 0; i < dimension; i++)
|
||||
v[i] = i < dim ? old[i] : 0;
|
||||
dim = dimension;
|
||||
delete [] old;
|
||||
}
|
||||
|
||||
void Vector::push_back(const double d) {
|
||||
resize(++dim);
|
||||
v[dim-1] = d;
|
||||
}
|
||||
|
||||
void Vector::set(const double d0, ...) {
|
||||
va_list args;
|
||||
v[0] = d0;
|
||||
va_start(args, d0);
|
||||
for (int i = 1; i < dim; i++)
|
||||
v[i] = va_arg(args, double);
|
||||
va_end(args);
|
||||
}
|
||||
|
||||
Vector& Vector::operator = (const Vector& dv) {
|
||||
if (this != &dv) {
|
||||
if (dim != dv.dim) {
|
||||
delete [] v;
|
||||
v = new double [dim = dv.dim];
|
||||
}
|
||||
for (int i = 0; i < dim; i++) v[i] = dv[i];
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
Vector& Vector::operator += (const Vector& dv) {
|
||||
for (int i = 0; i < dim; i++) v[i] += dv[i];
|
||||
return *this;
|
||||
}
|
||||
|
||||
Vector& Vector::operator -= (const Vector& dv) {
|
||||
for (int i = 0; i < dim; i++) v[i] -= dv[i];
|
||||
return *this;
|
||||
}
|
||||
|
||||
Vector& Vector::operator *= (double d) {
|
||||
for (int i = 0; i < dim; i++) v[i] *= d;
|
||||
return *this;
|
||||
}
|
||||
|
||||
Vector& Vector::operator /= (double d) {
|
||||
for (int i = 0; i < dim; i++) v[i] /= d;
|
||||
return *this;
|
||||
}
|
||||
|
||||
Vector operator - (const Vector& dv) {
|
||||
int dim = dv.dimension();
|
||||
Vector temp(dim);
|
||||
for (int i = 0; i < dim; i++) temp[i] = -dv[i];
|
||||
return temp;
|
||||
}
|
||||
|
||||
Vector operator * (const Vector& dv, double d) {
|
||||
int dim = dv.dimension();
|
||||
Vector temp(dim);
|
||||
for (int i = 0; i < dim; i++) temp[i] = dv[i] * d;
|
||||
return temp;
|
||||
}
|
||||
|
||||
Vector operator * (double d, const Vector& dv) {
|
||||
int dim = dv.dimension();
|
||||
Vector temp(dim);
|
||||
for (int i = 0; i < dim; i++) temp[i] = dv[i] * d;
|
||||
return temp;
|
||||
}
|
||||
|
||||
Vector operator / (const Vector& dv, double d) {
|
||||
int dim = dv.dimension();
|
||||
Vector temp(dim);
|
||||
for (int i = 0; i < dim; i++) temp[i] = dv[i] / d;
|
||||
return temp;
|
||||
}
|
||||
|
||||
Vector operator + (const Vector& v1, const Vector& v2) {
|
||||
int dim = v1.dimension();
|
||||
Vector temp(dim);
|
||||
for (int i = 0; i < dim; i++) temp[i] = v1[i] + v2[i];
|
||||
return temp;
|
||||
}
|
||||
|
||||
Vector operator - (const Vector& v1, const Vector& v2) {
|
||||
int dim = v1.dimension();
|
||||
Vector temp(dim);
|
||||
for (int i = 0; i < dim; i++) temp[i] = v1[i] - v2[i];
|
||||
return temp;
|
||||
}
|
||||
|
||||
double Vector::abs() {
|
||||
return std::sqrt(norm());
|
||||
}
|
||||
|
||||
double Vector::norm() {
|
||||
double sum = 0;
|
||||
for (int i = 0; i < dim; i++) sum += v[i] * v[i];
|
||||
return sum;
|
||||
}
|
||||
|
||||
double Vector::dot(const Vector& dv) {
|
||||
double sum = 0;
|
||||
for (int i = 0; i < dim; i++) sum += v[i] * dv[i];
|
||||
return sum;
|
||||
}
|
||||
|
||||
std::ostream& operator<<(std::ostream& os, const Vector& dv) {
|
||||
for (int i = 0; i < dv.dim; i++) {
|
||||
os << dv.v[i];
|
||||
if (i < dv.dim-1)
|
||||
os << '\t';
|
||||
else
|
||||
os << '\n';
|
||||
}
|
||||
return os;
|
||||
}
|
||||
|
||||
double dot(const Vector& v1, const Vector& v2) {
|
||||
double sum = 0;
|
||||
for (int i = 0; i < v1.size(); i++)
|
||||
sum += v1[i] * v2[i];
|
||||
return sum;
|
||||
}
|
||||
|
||||
} /* end namespace cpl */
|
||||
|
|
@ -0,0 +1,124 @@
|
|||
#ifndef CPL_VECTOR_HPP
|
||||
#define CPL_VECTOR_HPP
|
||||
|
||||
#include <iostream>
|
||||
|
||||
namespace cpl {
|
||||
|
||||
// Vector object with components of type double
|
||||
|
||||
class Vector {
|
||||
public:
|
||||
|
||||
Vector( // constructor, creates a new vector
|
||||
int dim = 1 // dimension, assumed = 1 if omitted
|
||||
); // components initialized to 0
|
||||
|
||||
Vector( // constructor, creates a new 2-d vector
|
||||
double v0, // value of first component
|
||||
double v1 // value of second componet
|
||||
);
|
||||
|
||||
Vector( // constructor, creates a new 3-d vector
|
||||
double v0, // value of first component
|
||||
double v1, // value of second component
|
||||
double v2 // value of third component
|
||||
);
|
||||
|
||||
Vector( // copy constructor, creates a new vector
|
||||
const Vector& // vector to be copied
|
||||
);
|
||||
|
||||
~Vector() { // destructor frees memory allocated by constructors
|
||||
delete [] v; // inline code to delete array v
|
||||
}
|
||||
|
||||
int dimension() const // returns number of components
|
||||
{ return dim; } // inline code to get dimension
|
||||
|
||||
int size() const // same as dimension()
|
||||
{ return dim; }
|
||||
|
||||
void resize( // resize the vector
|
||||
const int // new number of components
|
||||
); // old components preserved, any new zeroed
|
||||
|
||||
void push_back( // add a new component to the vector
|
||||
const double // value to set the new component
|
||||
);
|
||||
|
||||
void set( // set all the components of the vector
|
||||
const double, // value of first component
|
||||
... // second, third, ... , last component
|
||||
); // all components must be specified
|
||||
|
||||
const double // value returned by
|
||||
operator[]( // [ ] operator function
|
||||
const int i // index of component to return
|
||||
) const // value cannot be modified
|
||||
{ return v[i]; } // inline code to get value
|
||||
|
||||
double& // reference returned by
|
||||
operator[]( // [ ] operator function
|
||||
const int i // index of component to return
|
||||
) { return v[i]; } // inline code to get component
|
||||
|
||||
// operators to allow the vector to be assigned
|
||||
|
||||
Vector& operator = (const Vector&);
|
||||
|
||||
Vector& operator += (const Vector&);
|
||||
|
||||
Vector& operator -= (const Vector&);
|
||||
|
||||
// multiply or divide the vector by a scalar
|
||||
|
||||
Vector& operator *= (double);
|
||||
|
||||
Vector& operator /= (double);
|
||||
|
||||
double abs(); // return the length of the vector
|
||||
|
||||
double norm(); // return the norm of the vector
|
||||
|
||||
double dot( // return dot product of the vector
|
||||
const Vector& // with this vector
|
||||
);
|
||||
|
||||
// allow the << operator to be used to output vector components
|
||||
friend std::ostream& operator<<(std::ostream& os, const Vector& dv);
|
||||
|
||||
private:
|
||||
int dim; // store the number of components of the vector
|
||||
double *v; // pointer to stored components values
|
||||
};
|
||||
|
||||
// various arithmetic operations on vectors
|
||||
|
||||
// unary plus using inline code, just returns the same vector
|
||||
inline Vector operator + (const Vector& vec) { return vec; }
|
||||
|
||||
// unary minus, flips the sign of each component
|
||||
extern Vector operator - (const Vector&);
|
||||
|
||||
// multiplying and dividing a vector by a scalar
|
||||
|
||||
extern Vector operator * (const Vector&v, double);
|
||||
|
||||
extern Vector operator * (double, const Vector&);
|
||||
|
||||
extern Vector operator / (const Vector&, double);
|
||||
|
||||
// adding and subtracting two vectors
|
||||
|
||||
extern Vector operator + (const Vector&, const Vector&);
|
||||
|
||||
extern Vector operator - (const Vector&, const Vector&);
|
||||
|
||||
// dot product of two vectors
|
||||
|
||||
extern double dot(const Vector&, const Vector&);
|
||||
|
||||
} /* end namespace cpl */
|
||||
|
||||
#endif /* CPL_VECTOR_HPP */
|
||||
|
|
@ -0,0 +1,30 @@
|
|||
PROG = sho
|
||||
|
||||
SRCS = sho.cpp
|
||||
OBJS = ${SRCS:.cpp=.o}
|
||||
HDRS =
|
||||
|
||||
DEPS = popt
|
||||
CHEZ = /usr/lib/csv10.0.0/ta6le
|
||||
CPPFLAGS += -D_FORTIFY_SOURCE=3 -DCHEZDIR=\"${CHEZ}\"
|
||||
CXXFLAGS += -std=c++11 -Wall -Wextra -I${CHEZ} `pkg-config --cflags ${DEPS}`
|
||||
LDLIBS = ${LDFLAGS} -lm `pkg-config --libs ${DEPS}`
|
||||
|
||||
all: ${PROG}
|
||||
|
||||
${OBJS}: Makefile
|
||||
|
||||
${PROG}: ${OBJS}
|
||||
${CXX} -o $@ ${OBJS} ${LDLIBS} # ${CHEZ}/kernel.o
|
||||
|
||||
${SRCS}: ${HDRS}
|
||||
|
||||
lint:
|
||||
clang-tidy ${SRCS} ${HDRS} -- ${CPPFLAGS} ${CXXFLAGS}
|
||||
|
||||
plot: plot1.png
|
||||
|
||||
plot1.png: test.txt plot1.gnuplot
|
||||
gnuplot plot1.gnuplot
|
||||
|
||||
.PHONY: all lint
|
||||
|
|
@ -0,0 +1,46 @@
|
|||
#include <scheme.h>
|
||||
|
||||
#include "actkbd.h"
|
||||
|
||||
/* is the device currently grabbed? */
|
||||
static int grabbed_p(void) {
|
||||
/* this value comes from C, so an extra assert won't hurt */
|
||||
assert((grabbed == 0) || (grabbed == 1));
|
||||
return grabbed;
|
||||
}
|
||||
|
||||
void extension_init(void) {
|
||||
FILE *stream;
|
||||
|
||||
assert(config != NULL);
|
||||
|
||||
Sscheme_init(0);
|
||||
Sregister_boot_file(CHEZDIR "/petite.boot");
|
||||
Sregister_boot_file(CHEZDIR "/scheme.boot");
|
||||
Sbuild_heap(0, 0);
|
||||
|
||||
Sforeign_symbol("bind", (ptr) &set_cb);
|
||||
Sforeign_symbol("grab", (ptr) &set_grab);
|
||||
Sforeign_symbol("grabbed?", (ptr) &grabbed_p);
|
||||
Sforeign_symbol("led", (ptr) &set_led);
|
||||
|
||||
/* test whether we can read the config */
|
||||
stream = fopen(config, "r");
|
||||
if (!stream) {
|
||||
fprintf(stderr, "ERROR : %s not found or reading is not allowed.\n",
|
||||
config);
|
||||
exit(CONFERR);
|
||||
}
|
||||
fclose(stream);
|
||||
Sscheme_program(config, 0, NULL);
|
||||
}
|
||||
|
||||
/* a shortcut to reload stuff */
|
||||
void extension_reload(void) {
|
||||
Sscheme_program(config, 0, NULL);
|
||||
}
|
||||
|
||||
void extension_cleanup_callback(CB fun) {
|
||||
/* so the gc will be able to collect it */
|
||||
Sunlock_object(fun);
|
||||
}
|
||||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
|
|
@ -0,0 +1,19 @@
|
|||
set encoding utf8
|
||||
set terminal png size 800,600 enhanced font 'IBM Plex Sans Text,14'
|
||||
set output "plot1.png"
|
||||
|
||||
set xlabel 'idő [ms]'
|
||||
set ylabel 'kitérés [m]'
|
||||
|
||||
|
||||
set style line 100 lt 1 lc rgb "grey" lw 0.5
|
||||
set grid ls 100
|
||||
|
||||
set style line 103 lw 4 lt rgb "#b8ee30"
|
||||
set style line 104 pt 1 ps 1 lw 1 lt rgb "#8b1c62"
|
||||
set title 'Harmonikus oszcillátor'
|
||||
|
||||
plot 'test.txt' title 'kitérés' with points ls 104
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,77 @@
|
|||
#include <cmath>
|
||||
#include <iostream>
|
||||
#include <functional>
|
||||
|
||||
#include <popt.h>
|
||||
|
||||
using namespace std;
|
||||
|
||||
double omega; // the natural frequency
|
||||
double x, v; // position and velocity at time t
|
||||
int periods; // number of periods to integrate
|
||||
int steps; // number of time steps dt per period
|
||||
|
||||
void getInput(void); // for user to input parameters
|
||||
void EulerCromer(double dt); // takes an Euler-Cromer step
|
||||
void Euler(double dt); // takes an Euler step
|
||||
double energy(void); // computes the energy
|
||||
void sim(function<void(double)> integrator); // runs the simulation
|
||||
|
||||
typedef double v2d __attribute__ ((vector_size (2*8)));
|
||||
|
||||
enum {
|
||||
OK=0,
|
||||
USAGE=1,
|
||||
HOSTFAIL=2,
|
||||
};
|
||||
|
||||
void EulerCromer (double dt) {
|
||||
double a = - omega * omega * x;
|
||||
v += a * dt;
|
||||
x += v * dt;
|
||||
}
|
||||
|
||||
void Euler(double dt) {
|
||||
double a = - omega * omega * x;
|
||||
x += v * dt;
|
||||
v += a * dt;
|
||||
}
|
||||
|
||||
double energy(void) {
|
||||
return 0.5 * (v * v + omega * omega * x * x);
|
||||
}
|
||||
|
||||
void sim(function<void(double)> integrator) {
|
||||
const double pi = 4 * atan(1.0);
|
||||
double T = 2 * pi / omega;
|
||||
double dt = T / steps;
|
||||
double t = 0;
|
||||
cout << t << '\t' << x << '\t' << v << '\t' << energy() << '\n';
|
||||
for (int p = 1; p <= periods; p++) {
|
||||
for (int s = 0; s < steps; s++) {
|
||||
integrator(dt);
|
||||
t += dt;
|
||||
cout << t << '\t' << x << '\t' << v << '\t' << energy() << '\n';
|
||||
}
|
||||
// cout << "Period = " << p << "\tt = " << t
|
||||
// cout << "\tx = " << x << "\tv = " << v
|
||||
// cout << "\tenergy = " << energy() << endl;
|
||||
}
|
||||
}
|
||||
|
||||
int main (int argc, const char *argv[]) {
|
||||
const struct poptOption options_table[] = {
|
||||
{"omega", 'o', POPT_ARG_DOUBLE, &omega, 0, "Value of omega", NULL},
|
||||
{"xnull", 'x', POPT_ARG_DOUBLE, &x, 0, "Value of x(0)", NULL},
|
||||
{"vnull", 'v', POPT_ARG_DOUBLE, &v, 0, "Value of v(0)", NULL},
|
||||
{"periods", 'p', POPT_ARG_INT, &periods, 0, "Number of periods", NULL},
|
||||
{"steps", 's', POPT_ARG_INT, &steps, 0, "Steps per period", NULL},
|
||||
POPT_AUTOHELP
|
||||
{NULL, 0, 0, NULL, 0, NULL, NULL}};
|
||||
|
||||
poptContext opt_con = poptGetContext("sho", argc, argv, options_table, 0);
|
||||
poptSetOtherOptionHelp(opt_con, "[OPTIONS]*");
|
||||
poptGetNextOpt(opt_con);
|
||||
sim(Euler);
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
Binary file not shown.
|
|
@ -0,0 +1,201 @@
|
|||
0 1 0.1 50.005
|
||||
0.00628319 0.99668 -0.528319 49.8082
|
||||
0.0125664 0.989426 -1.15455 49.6147
|
||||
0.0188496 0.978266 -1.77623 49.4277
|
||||
0.0251327 0.963243 -2.39089 49.2501
|
||||
0.0314159 0.944418 -2.99611 49.0846
|
||||
0.0376991 0.921865 -3.58951 48.934
|
||||
0.0439823 0.895672 -4.16873 48.8006
|
||||
0.0502655 0.865943 -4.7315 48.6864
|
||||
0.0565487 0.832795 -5.27559 48.5933
|
||||
0.0628319 0.79636 -5.79885 48.5228
|
||||
0.069115 0.756781 -6.29922 48.4759
|
||||
0.0753982 0.714214 -6.77472 48.4535
|
||||
0.0816814 0.668828 -7.22347 48.4558
|
||||
0.0879646 0.620801 -7.64371 48.4828
|
||||
0.0942478 0.570323 -8.03377 48.5342
|
||||
0.100531 0.517594 -8.39211 48.609
|
||||
0.106814 0.462822 -8.71733 48.7061
|
||||
0.113097 0.406222 -9.00813 48.824
|
||||
0.119381 0.348018 -9.26336 48.9608
|
||||
0.125664 0.288441 -9.48203 49.1144
|
||||
0.131947 0.227725 -9.66326 49.2823
|
||||
0.13823 0.16611 -9.80635 49.4618
|
||||
0.144513 0.103839 -9.91072 49.6503
|
||||
0.150796 0.0411583 -9.97596 49.8446
|
||||
0.15708 -0.021685 -10.0018 50.0417
|
||||
0.163363 -0.0844427 -9.9882 50.2386
|
||||
0.169646 -0.146867 -9.93514 50.432
|
||||
0.175929 -0.208712 -9.84286 50.619
|
||||
0.182212 -0.269732 -9.71172 50.7965
|
||||
0.188496 -0.329688 -9.54224 50.9619
|
||||
0.194779 -0.388342 -9.3351 51.1125
|
||||
0.201062 -0.445463 -9.09109 51.2458
|
||||
0.207345 -0.500825 -8.8112 51.3599
|
||||
0.213628 -0.554211 -8.49652 51.4529
|
||||
0.219911 -0.605408 -8.1483 51.5233
|
||||
0.226195 -0.654215 -7.76791 51.5701
|
||||
0.232478 -0.70044 -7.35686 51.5925
|
||||
0.238761 -0.743899 -6.91676 51.59
|
||||
0.245044 -0.784421 -6.44935 51.5629
|
||||
0.251327 -0.821847 -5.95649 51.5115
|
||||
0.257611 -0.856028 -5.4401 51.4366
|
||||
0.263894 -0.88683 -4.90225 51.3394
|
||||
0.270177 -0.914131 -4.34503 51.2214
|
||||
0.27646 -0.937822 -3.77067 51.0845
|
||||
0.282743 -0.957812 -3.18142 50.9309
|
||||
0.289027 -0.97402 -2.57961 50.7629
|
||||
0.29531 -0.986383 -1.96761 50.5833
|
||||
0.301593 -0.994852 -1.34785 50.3948
|
||||
0.307876 -0.999393 -0.722766 50.2005
|
||||
0.314159 -0.999989 -0.094829 50.0034
|
||||
0.320442 -0.996637 0.533482 49.8065
|
||||
0.326726 -0.98935 1.15969 49.6131
|
||||
0.333009 -0.978158 1.78131 49.4262
|
||||
0.339292 -0.963104 2.39591 49.2487
|
||||
0.345575 -0.944248 3.00105 49.0833
|
||||
0.351858 -0.921664 3.59433 48.9328
|
||||
0.358142 -0.895442 4.17343 48.7996
|
||||
0.364425 -0.865684 4.73606 48.6856
|
||||
0.370708 -0.832509 5.27998 48.5927
|
||||
0.376991 -0.796047 5.80306 48.5223
|
||||
0.383274 -0.756443 6.30323 48.4757
|
||||
0.389557 -0.713852 6.77852 48.4534
|
||||
0.395841 -0.668443 7.22705 48.4559
|
||||
0.402124 -0.620396 7.64704 48.4832
|
||||
0.408407 -0.569899 8.03685 48.5347
|
||||
0.41469 -0.517152 8.39493 48.6097
|
||||
0.420973 -0.462363 8.71986 48.707
|
||||
0.427257 -0.405749 9.01037 48.825
|
||||
0.43354 -0.347534 9.26531 48.962
|
||||
0.439823 -0.287946 9.48367 49.1157
|
||||
0.446106 -0.227221 9.6646 49.2837
|
||||
0.452389 -0.1656 9.80736 49.4634
|
||||
0.458673 -0.103325 9.91141 49.6519
|
||||
0.464956 -0.0406416 9.97633 49.8462
|
||||
0.471239 0.022202 10.0019 50.0433
|
||||
0.477522 0.0849579 9.98792 50.2402
|
||||
0.483805 0.147379 9.93454 50.4336
|
||||
0.490088 0.209217 9.84194 50.6205
|
||||
0.496372 0.27023 9.71048 50.798
|
||||
0.502655 0.330176 9.54069 50.9632
|
||||
0.508938 0.388818 9.33324 51.1136
|
||||
0.515221 0.445926 9.08894 51.2469
|
||||
0.521504 0.501273 8.80875 51.3608
|
||||
0.527788 0.554641 8.49379 51.4536
|
||||
0.534071 0.605819 8.1453 51.5238
|
||||
0.540354 0.654606 7.76465 51.5704
|
||||
0.546637 0.700809 7.35335 51.5925
|
||||
0.55292 0.744244 6.91302 51.5899
|
||||
0.559203 0.784742 6.4454 51.5626
|
||||
0.565487 0.822142 5.95233 51.511
|
||||
0.57177 0.856296 5.43576 51.4359
|
||||
0.578053 0.887069 4.89774 51.3385
|
||||
0.584336 0.91434 4.34038 51.2204
|
||||
0.590619 0.938002 3.76588 51.0833
|
||||
0.596903 0.957961 3.17652 50.9296
|
||||
0.603186 0.974137 2.57461 50.7615
|
||||
0.609469 0.986469 1.96254 50.5818
|
||||
0.615752 0.994905 1.34273 50.3933
|
||||
0.622035 0.999414 0.717608 50.1989
|
||||
0.628319 0.999977 0.0896581 50.0018
|
||||
0.634602 0.996593 -0.538646 49.8049
|
||||
0.640885 0.989274 -1.16482 49.6116
|
||||
0.647168 0.97805 -1.7864 49.4247
|
||||
0.653451 0.962964 -2.40093 49.2472
|
||||
0.659734 0.944077 -3.00598 49.082
|
||||
0.666018 0.921463 -3.59916 48.9317
|
||||
0.672301 0.895211 -4.17813 48.7985
|
||||
0.678584 0.865425 -4.74061 48.6847
|
||||
0.684867 0.832222 -5.28437 48.592
|
||||
0.69115 0.795734 -5.80727 48.5218
|
||||
0.697434 0.756104 -6.30725 48.4754
|
||||
0.703717 0.71349 -6.78232 48.4533
|
||||
0.71 0.668059 -7.23062 48.456
|
||||
0.716283 0.61999 -7.65037 48.4835
|
||||
0.722566 0.569473 -8.03992 48.5352
|
||||
0.728849 0.516709 -8.39774 48.6104
|
||||
0.735133 0.461905 -8.72239 48.7079
|
||||
0.741416 0.405277 -9.01262 48.8261
|
||||
0.747699 0.347049 -9.26726 48.9632
|
||||
0.753982 0.287451 -9.48532 49.117
|
||||
0.760265 0.226718 -9.66593 49.2851
|
||||
0.766549 0.16509 -9.80838 49.4649
|
||||
0.772832 0.10281 -9.91211 49.6534
|
||||
0.779115 0.0401249 -9.9767 49.8478
|
||||
0.785398 -0.022719 -10.0019 50.045
|
||||
0.791681 -0.0854732 -9.98764 50.2418
|
||||
0.797965 -0.14789 -9.93394 50.4351
|
||||
0.804248 -0.209723 -9.84102 50.622
|
||||
0.810531 -0.270728 -9.70924 50.7994
|
||||
0.816814 -0.330664 -9.53914 50.9645
|
||||
0.823097 -0.389295 -9.33138 51.1148
|
||||
0.82938 -0.446389 -9.08678 51.2479
|
||||
0.835664 -0.50172 -8.8063 51.3616
|
||||
0.841947 -0.555071 -8.49106 51.4543
|
||||
0.84823 -0.606231 -8.1423 51.5243
|
||||
0.854513 -0.654997 -7.76139 51.5707
|
||||
0.860796 -0.701178 -7.34985 51.5926
|
||||
0.86708 -0.74459 -6.90928 51.5898
|
||||
0.873363 -0.785063 -6.44144 51.5623
|
||||
0.879646 -0.822436 -5.94817 51.5104
|
||||
0.885929 -0.856563 -5.43142 51.4352
|
||||
0.892212 -0.887308 -4.89323 51.3376
|
||||
0.898495 -0.91455 -4.33572 51.2193
|
||||
0.904779 -0.938182 -3.76109 51.0821
|
||||
0.911062 -0.958109 -3.17161 50.9282
|
||||
0.917345 -0.974255 -2.56961 50.7601
|
||||
0.923628 -0.986554 -1.95747 50.5803
|
||||
0.929911 -0.994958 -1.3376 50.3917
|
||||
0.936195 -0.999435 -0.71245 50.1973
|
||||
0.942478 -0.999966 -0.0844871 50.0001
|
||||
0.948761 -0.996549 0.54381 49.8033
|
||||
0.955044 -0.989198 1.16996 49.61
|
||||
0.961327 -0.977941 1.79149 49.4232
|
||||
0.967611 -0.962824 2.40595 49.2458
|
||||
0.973894 -0.943906 3.01091 49.0807
|
||||
0.980177 -0.921262 3.60398 48.9305
|
||||
0.98646 -0.89498 4.18283 48.7975
|
||||
0.992743 -0.865166 4.74516 48.6839
|
||||
0.999026 -0.831935 5.28876 48.5913
|
||||
1.00531 -0.795421 5.81148 48.5214
|
||||
1.01159 -0.755766 6.31126 48.4751
|
||||
1.01788 -0.713127 6.78612 48.4533
|
||||
1.02416 -0.667674 7.23419 48.4562
|
||||
1.03044 -0.619584 7.6537 48.4838
|
||||
1.03673 -0.569048 8.043 48.5357
|
||||
1.04301 -0.516266 8.40054 48.6111
|
||||
1.04929 -0.461446 8.72492 48.7088
|
||||
1.05558 -0.404804 9.01486 48.8271
|
||||
1.06186 -0.346564 9.2692 48.9644
|
||||
1.06814 -0.286955 9.48696 49.1183
|
||||
1.07442 -0.226214 9.66726 49.2866
|
||||
1.08071 -0.16458 9.80939 49.4664
|
||||
1.08699 -0.102296 9.9128 49.655
|
||||
1.09327 -0.0396082 9.97707 49.8494
|
||||
1.09956 0.023236 10.002 50.0466
|
||||
1.10584 0.0859884 9.98736 50.2434
|
||||
1.11212 0.148401 9.93333 50.4367
|
||||
1.11841 0.210228 9.84009 50.6235
|
||||
1.12469 0.271226 9.708 50.8008
|
||||
1.13097 0.331152 9.53758 50.9658
|
||||
1.13726 0.389771 9.32951 51.116
|
||||
1.14354 0.446851 9.08461 51.2489
|
||||
1.14982 0.502168 8.80385 51.3625
|
||||
1.15611 0.555501 8.48833 51.4549
|
||||
1.16239 0.606642 8.13929 51.5248
|
||||
1.16867 0.655388 7.75813 51.571
|
||||
1.17496 0.701546 7.34634 51.5927
|
||||
1.18124 0.744935 6.90554 51.5897
|
||||
1.18752 0.785383 6.43749 51.5619
|
||||
1.19381 0.82273 5.94402 51.5099
|
||||
1.20009 0.85683 5.42708 51.4345
|
||||
1.20637 0.887546 4.88872 51.3367
|
||||
1.21265 0.914759 4.33106 51.2182
|
||||
1.21894 0.938361 3.7563 51.0809
|
||||
1.22522 0.958258 3.16671 50.9269
|
||||
1.2315 0.974372 2.56462 50.7586
|
||||
1.23779 0.986639 1.9524 50.5788
|
||||
1.24407 0.995011 1.33248 50.3901
|
||||
1.25035 0.999455 0.707292 50.1957
|
||||
1.25664 0.999954 0.0793161 49.9985
|
||||
|
|
@ -0,0 +1,30 @@
|
|||
PROG = pen
|
||||
|
||||
SRCS = pendulum.cpp pendulum-gl.cpp
|
||||
OBJS = ${SRCS:.cpp=.o}
|
||||
HDRS =
|
||||
|
||||
DEPS = popt gl glu glut
|
||||
CHEZ = /usr/lib/csv10.0.0/ta6le
|
||||
CPPFLAGS += -D_FORTIFY_SOURCE=3 -DCHEZDIR=\"${CHEZ}\"
|
||||
CXXFLAGS += -std=c++11 -Wall -Wextra -I${CHEZ} `pkg-config --cflags ${DEPS}`
|
||||
LDLIBS = ${LDFLAGS} -lm `pkg-config --libs ${DEPS}`
|
||||
|
||||
all: ${PROG}
|
||||
|
||||
${OBJS}: Makefile
|
||||
|
||||
${PROG}: ${OBJS}
|
||||
${CXX} -o $@ ${OBJS} ${LDLIBS} # ${CHEZ}/kernel.o
|
||||
|
||||
${SRCS}: ${HDRS}
|
||||
|
||||
lint:
|
||||
clang-tidy ${SRCS} ${HDRS} -- ${CPPFLAGS} ${CXXFLAGS}
|
||||
|
||||
plot: plot1.png
|
||||
|
||||
plot1.png: test.txt plot1.gnuplot
|
||||
gnuplot plot1.gnuplot
|
||||
|
||||
.PHONY: all lint
|
||||
|
|
@ -0,0 +1,10 @@
|
|||
#include "pendulum.hpp"
|
||||
|
||||
void adaptiveRK4Step(State x, double dt, double accuracy, State f);
|
||||
|
||||
|
||||
|
||||
void adaptiveRK4Step(State x, double dt, double accuracy, State f)
|
||||
{
|
||||
|
||||
}
|
||||
|
|
@ -0,0 +1,246 @@
|
|||
#include <cmath>
|
||||
#include <cstdlib>
|
||||
#include <iostream>
|
||||
#include <sstream>
|
||||
#include <string>
|
||||
using namespace std;
|
||||
|
||||
#ifdef __APPLE__ // Mac OS X uses different header
|
||||
# include <GLUT/glut.h>
|
||||
#else // Unix and Windows
|
||||
# include <GL/gl.h>
|
||||
# include <GL/glu.h>
|
||||
# include <GL/glut.h>
|
||||
#endif
|
||||
|
||||
#include "vector.hpp" // vectors with components of type double
|
||||
#include "odeint.hpp" // ODE integration routines, Runge-Kutta ...
|
||||
using namespace cpl;
|
||||
|
||||
const double pi = 4 * atan(1.0);
|
||||
|
||||
const double g = 9.8; // acceleration of gravity
|
||||
|
||||
double L = 1.0; // length of pendulum
|
||||
double q = 0.5; // damping coefficient
|
||||
double Omega_D = 2.0/3.0; // frequency of driving force
|
||||
double F_D = 0.9; // amplitude of driving force
|
||||
bool nonlinear; // linear if false
|
||||
|
||||
Vector f(const Vector& x) { // extended derivative vector
|
||||
double t = x[0];
|
||||
double theta = x[1];
|
||||
double omega = x[2];
|
||||
Vector f(3); // Vector with 3 components
|
||||
f[0] = 1;
|
||||
f[1] = omega;
|
||||
if (nonlinear)
|
||||
f[2] = - (g/L) * sin(theta) - q * omega + F_D * sin(Omega_D * t);
|
||||
else
|
||||
f[2] = - (g/L) * theta - q * omega + F_D * sin(Omega_D * t);
|
||||
return f;
|
||||
}
|
||||
|
||||
Vector x(3);
|
||||
|
||||
void getInput() {
|
||||
cout << " Nonlinear damped driven pendulum\n"
|
||||
<< " --------------------------------\n"
|
||||
<< " Enter linear or nonlinear: ";
|
||||
string response;
|
||||
cin >> response;
|
||||
nonlinear = (response[0] == 'n');
|
||||
cout<< " Length of pendulum L: ";
|
||||
cin >> L;
|
||||
cout<< " Enter damping coefficient q: ";
|
||||
cin >> q;
|
||||
cout << " Enter driving frequency Omega_D: ";
|
||||
cin >> Omega_D;
|
||||
cout << " Enter driving amplitude F_D: ";
|
||||
cin >> F_D;
|
||||
cout << " Enter theta(0) and omega(0): ";
|
||||
double theta, omega;
|
||||
cin >> theta >> omega;
|
||||
|
||||
x[0] = 0;
|
||||
x[1] = theta;
|
||||
x[2] = omega;
|
||||
}
|
||||
|
||||
double dt = 0.05;
|
||||
double accuracy = 1e-6;
|
||||
|
||||
void step() {
|
||||
adaptiveRK4Step(x, dt, accuracy, f);
|
||||
}
|
||||
|
||||
double frames_per_second = 30; // for animation in real time
|
||||
|
||||
void animation_step() {
|
||||
double start = x[0];
|
||||
clock_t start_time = clock();
|
||||
step();
|
||||
double tau = 1.0 / frames_per_second;
|
||||
while (x[0] - start < tau)
|
||||
step();
|
||||
while ((double(clock()) - start_time)/CLOCKS_PER_SEC < tau)
|
||||
;
|
||||
glutPostRedisplay();
|
||||
}
|
||||
|
||||
void drawText(const string& str, double x, double y) {
|
||||
glRasterPos2d(x, y);
|
||||
int len = str.find('\0');
|
||||
for (int i = 0; i < len; i++)
|
||||
glutBitmapCharacter(GLUT_BITMAP_HELVETICA_12, str[i]);
|
||||
}
|
||||
|
||||
void drawVariables() {
|
||||
// write t, theta, omega values
|
||||
glColor3ub(0, 0, 255);
|
||||
ostringstream os;
|
||||
os << "t = " << x[0] << ends;
|
||||
drawText(os.str(), -1.1, -1.1);
|
||||
os.seekp(0);
|
||||
os << "theta = " << x[1] << ends;
|
||||
drawText(os.str(), -1.1, 1.1);
|
||||
os.seekp(0);
|
||||
os << "omega = " << x[2] << ends;
|
||||
drawText(os.str(), 0.3, 1.1);
|
||||
}
|
||||
|
||||
void displayPendulum() {
|
||||
glClear(GL_COLOR_BUFFER_BIT);
|
||||
|
||||
// draw the pendulum rod
|
||||
double xp = sin(x[1]);
|
||||
double yp = -cos(x[1]);
|
||||
glColor3ub(0, 255, 0);
|
||||
glBegin(GL_LINES);
|
||||
glVertex2d(0, 0);
|
||||
glVertex2d(xp, yp);
|
||||
glEnd();
|
||||
|
||||
// draw the pendulum bob
|
||||
glPushMatrix();
|
||||
glTranslated(sin(x[1]), -cos(x[1]), 0);
|
||||
glColor3ub(255, 0, 0);
|
||||
const double r = 0.1;
|
||||
glPolygonMode(GL_FRONT, GL_FILL);
|
||||
glBegin(GL_POLYGON);
|
||||
for (int i = 0; i < 12; i++) {
|
||||
double theta = 2 * pi * i / 12.0;
|
||||
glVertex2d(r * cos(theta), r * sin(theta));
|
||||
}
|
||||
glEnd();
|
||||
glPopMatrix();
|
||||
|
||||
// write t, theta, and omega
|
||||
drawVariables();
|
||||
|
||||
// we are using double buffering - write buffer to screen
|
||||
glutSwapBuffers();
|
||||
}
|
||||
|
||||
bool clearPhasePlot;
|
||||
|
||||
void displayPhasePlot() {
|
||||
static double thetaOld, omegaOld;
|
||||
if (clearPhasePlot) {
|
||||
glClear(GL_COLOR_BUFFER_BIT);
|
||||
clearPhasePlot = false;
|
||||
thetaOld = 2 * pi, omegaOld = 2 * pi;
|
||||
|
||||
// draw axes
|
||||
glColor3ub(0, 255, 0);
|
||||
glBegin(GL_LINES);
|
||||
glVertex2d(0, -1); glVertex2d(0, 1);
|
||||
glVertex2d(-1, 0); glVertex2d(1, 0);
|
||||
glEnd();
|
||||
}
|
||||
|
||||
// draw the phase point
|
||||
double theta = x[1];
|
||||
while (theta >= pi) theta -= 2 * pi;
|
||||
while (theta < -pi) theta += 2 * pi;
|
||||
double omega = x[2];
|
||||
glColor3ub(255, 0, 0);
|
||||
if (abs(theta - thetaOld) < pi && abs(omega) < pi
|
||||
&& abs(omega - omegaOld) < pi) {
|
||||
glBegin(GL_LINES);
|
||||
glVertex2d(thetaOld / pi, omegaOld / pi);
|
||||
glVertex2d(theta / pi, omega / pi);
|
||||
glEnd();
|
||||
}
|
||||
thetaOld = theta, omegaOld = omega;
|
||||
glPopMatrix();
|
||||
|
||||
// write t, theta, and omega
|
||||
glColor3ub(255, 255, 255);
|
||||
glRectd(-1.2, 1, 1.2, 1.2);
|
||||
glRectd(-1.2, -1, 1.2, -1.2);
|
||||
drawVariables();
|
||||
|
||||
glutSwapBuffers();
|
||||
}
|
||||
|
||||
void reshape(int w, int h) {
|
||||
glViewport(0, 0, w, h);
|
||||
glMatrixMode(GL_PROJECTION);
|
||||
glLoadIdentity();
|
||||
double aspectRatio = w/double(h);
|
||||
double d = 1.2;
|
||||
if (aspectRatio > 1)
|
||||
glOrtho(-d*aspectRatio, d*aspectRatio, -d, d, -1.0, 1.0);
|
||||
else
|
||||
glOrtho(-d, d, -d/aspectRatio, d/aspectRatio, -1.0, 1.0);
|
||||
glMatrixMode(GL_MODELVIEW);
|
||||
}
|
||||
|
||||
bool running; // is the animation running?
|
||||
bool phasePlot; // switch to phase plot if true
|
||||
|
||||
void mouse(int button, int state, int x, int y) {
|
||||
switch (button) {
|
||||
case GLUT_LEFT_BUTTON:
|
||||
if (state == GLUT_DOWN) {
|
||||
if (running) {
|
||||
glutIdleFunc(NULL);
|
||||
running = false;
|
||||
} else {
|
||||
glutIdleFunc(animation_step);
|
||||
running = true;
|
||||
}
|
||||
}
|
||||
break;
|
||||
case GLUT_RIGHT_BUTTON:
|
||||
if (state == GLUT_DOWN) {
|
||||
if (phasePlot) {
|
||||
glutDisplayFunc(displayPendulum);
|
||||
phasePlot = false;
|
||||
} else {
|
||||
glutDisplayFunc(displayPhasePlot);
|
||||
clearPhasePlot = phasePlot = true;
|
||||
}
|
||||
glutPostRedisplay();
|
||||
}
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
int main(int argc, char *argv[]) {
|
||||
glutInit(&argc, argv);
|
||||
getInput();
|
||||
glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB);
|
||||
glutInitWindowSize(600, 600);
|
||||
glutInitWindowPosition(100, 100);
|
||||
glutCreateWindow(argv[0]);
|
||||
glClearColor(1.0, 1.0, 1.0, 0.0);
|
||||
glShadeModel(GL_FLAT);
|
||||
glutDisplayFunc(displayPendulum);
|
||||
glutReshapeFunc(reshape);
|
||||
glutMouseFunc(mouse);
|
||||
glutIdleFunc(NULL);
|
||||
glutMainLoop();
|
||||
}
|
||||
|
||||
|
|
@ -0,0 +1,76 @@
|
|||
#include <cmath>
|
||||
#include <cstdlib>
|
||||
#include <iostream>
|
||||
|
||||
#include <popt.h>
|
||||
|
||||
using namespace std;
|
||||
|
||||
#include "pendulum.hpp"
|
||||
#include "odeint.hpp" // ODE integration routines, Runge-Kutta ...
|
||||
|
||||
using namespace cpl;
|
||||
|
||||
const double pi = 4 * atan(1.0);
|
||||
|
||||
const double g = 9.8; // acceleration of gravity
|
||||
|
||||
double L = 1.0; // length of pendulum
|
||||
double q = 0.5; // damping coefficient
|
||||
double Omega_D = 2.0/3.0; // frequency of driving force
|
||||
double F_D = 0.9; // amplitude of driving force
|
||||
bool nonlinear; // linear if false
|
||||
|
||||
State f(const State& s) { // extended derivative vector
|
||||
double t = s.t;
|
||||
double theta = s.theta;
|
||||
double omega = s.omega;
|
||||
State f; // Vector with 3 components
|
||||
f.t = 1;
|
||||
|
||||
f.theta = omega;
|
||||
|
||||
f.omega = (- (g/L) * (nonlinear ? sin(theta) : theta ) -
|
||||
q * omega + F_D * sin(Omega_D * t));
|
||||
|
||||
return f;
|
||||
}
|
||||
|
||||
int main() {
|
||||
const struct poptOption options_table[] = {
|
||||
{"length", 'l', POPT_ARG_DOUBLE, &L, 0, "Length of pendulum", NULL},
|
||||
{"damping", 'q', POPT_ARG_DOUBLE, &q, 0, "Damping coefficient", NULL},
|
||||
{"omegaD", 'o', POPT_ARG_DOUBLE, &Omega_D, 0, "Value of Omega_D", NULL},
|
||||
{"F_D", 'f', POPT_ARG_DOUBLE, &F_D, 0, "Value of F_D", NULL},
|
||||
{"theta_0", 't', POPT_ARG_DOUBLE, &theta, 0, "Value of theta(0)", NULL},
|
||||
{"omega_0", '0', POPT_ARG_DOUBLE, &omega, 0, "Value of omega(0)", NULL},
|
||||
{"tmax", 'm', POPT_ARG_DOUBLE, &tMax, 0, "Integration time", NULL},
|
||||
POPT_AUTOHELP
|
||||
{NULL, 0, 0, NULL, 0, NULL, NULL}};
|
||||
|
||||
poptContext opt_con = poptGetContext("sho", argc, argv, options_table, 0);
|
||||
poptSetOtherOptionHelp(opt_con, "[OPTIONS]*");
|
||||
poptGetNextOpt(opt_con);
|
||||
|
||||
double dt = 0.05;
|
||||
double accuracy = 1e-6;
|
||||
|
||||
double t = 0;
|
||||
State x;
|
||||
x.t = t;
|
||||
x.theta = theta;
|
||||
x.omega = omega;
|
||||
|
||||
while (t < tMax) {
|
||||
adaptiveRK4Step(x, dt, accuracy, f);
|
||||
t = x.t, theta = x.theta, omega = x;
|
||||
if (nonlinear) {
|
||||
while (theta >= pi) theta -= 2 * pi;
|
||||
while (theta < -pi) theta += 2 * pi;
|
||||
}
|
||||
cout << t << '\t' << theta << '\t' << omega << '\n';
|
||||
}
|
||||
|
||||
cerr << "Finished." << endl;
|
||||
}
|
||||
|
||||
|
|
@ -0,0 +1,9 @@
|
|||
#ifndef PENDULUM_H
|
||||
#define PENDULUM_H
|
||||
struct State {
|
||||
double t;
|
||||
double theta;
|
||||
double omega;
|
||||
};
|
||||
|
||||
#endif
|
||||
|
|
@ -0,0 +1,19 @@
|
|||
set encoding utf8
|
||||
set terminal png size 800,600 enhanced font 'IBM Plex Sans Text,14'
|
||||
set output "plot1.png"
|
||||
|
||||
set xlabel 'idő [ms]'
|
||||
set ylabel 'kitérés [m]'
|
||||
|
||||
|
||||
set style line 100 lt 1 lc rgb "grey" lw 0.5
|
||||
set grid ls 100
|
||||
|
||||
set style line 103 lw 4 lt rgb "#b8ee30"
|
||||
set style line 104 pt 1 ps 1 lw 1 lt rgb "#8b1c62"
|
||||
set title 'Harmonikus oszcillátor'
|
||||
|
||||
plot 'test.txt' title 'kitérés' with points ls 104
|
||||
|
||||
|
||||
|
||||
Loading…
Reference in New Issue