---

_______________________＿＿＿＿__＿＿＿_____＿＿＿________計算コードの具体例_______________________________________________
1、三次方程式(コード1)
#include <iostream>
#include <string>
#include <cmath>
#include <complex>
#include <vector>
#include <fstream>
#include <algorithm>//transform
using namespace std;

const double pi=3.14159265358979;
const complex<double> real_unit(1.0,0.0),imag_unit(0.0,1.0);

complex<double> csum0(complex<double> d1,complex<double> d2){;return (pow(d1,1/3.0)+pow(d2,1/3.0));}

class binary_csum1
{
public: complex<double> operator()(complex<double> x,complex<double> y)
{
return pow(x,1/3.0)+pow(y,1/3.0);
}
};

int main(void)
{
double x1=12.0,x2=2.0,x3=8.0;//解(真値)
ofstream out0;
out0.open("in3eq.csv", ios::out | ios::binary );
out0<<-x1*x2*x3<<endl;
out0<<(x1*x2+x1*x3+x2*x3)<<endl;
out0<<-(x1+x2+x3)<<endl;
out0.close();

string string_data0;//データを文字列で読み、strtodで数値に変換し、vectorに格納する。
vector<double>a;//a[0],a[1],a[2]となる。x の次数の係数である。
double indata;

ifstream input0("in3eq.csv", ios::in | ios::binary);//
if(!input0.is_open())
{
cout<<"Cannot open file　\n";exit(1);
}
while(!input0.eof()) //eof()関数で最後までデータを読込む。
{
input0>>string_data0;
if(!input0.eof())
{
indata=strtod((const char*)string_data0.c_str(),NULL);//文字列を数値に変更。
a.push_back(indata);
}
}
input0.close();
a.resize(a.size());

//x^3+a[2]*x^2+a[1]*x+a[0];
double p,q;
p=a[1]-a[2]*a[2]/3.0;
q=a[0]-a[1]*a[2]/3.0+2.0*a[2]*a[2]*a[2]/27.0;

complex<double>z(0.0,0.0);
z=(q*q/4.0+p*p*p/27.0)*real_unit;
complex<double> z1(0.0,0.0);
z1=pow(z,1/2.0);
complex<double> y1(0.0,0.0),y2(0.0,0.0);
y1=-q/2.0+z1;
y2=-q/2.0-z1;

complex<double>c0(0.0,0.0);
transform(&y1,&y1+1,&y2,&c0,csum0);//transform使用例。

complex<double>c1(0.0,0.0);
transform(&y1,&y1+1,&y2,&c1,binary_csum1());//transform使用例。

ofstream outk0;
outk0.open("data_check_with_transform.csv", ios::out | ios::binary );
outk0<<"complex y1= "<<","<<y1<<endl;
outk0<<"complex y2= "<<","<<y2<<endl;
outk0<<" c0 of transform csum1 = y1^1/3+y2^1/3 = "<<","<<c0<<endl;
outk0<<" c1 of transform binary_csum2() = y1^1/3+y2^1/3 = "<<","<<c1<<endl;
outk0.close();


complex<double> shift(-1.0/2.0,sqrt(3.0)/2.0);
vector<complex<double>>y(9,(0.0,0.0));//複素数解を含めると9個、必要になる。

y[0]=1.0*pow(y1,1/3.0)+1.0*pow(y2,1/3.0);
y[1]=1.0*pow(y1,1/3.0)+shift*pow(y2,1/3.0);
y[2]=1.0*pow(y1,1/3.0)+shift*shift*pow(y2,1/3.0);

y[3]=shift*pow(y1,1/3.0)+1.0*pow(y2,1/3.0);
y[4]=shift*pow(y1,1/3.0)+shift*pow(y2,1/3.0);
y[5]=shift*pow(y1,1/3.0)+shift*shift*pow(y2,1/3.0);

y[6]=shift*shift*pow(y1,1/3.0)+1.0*pow(y2,1/3.0);
y[7]=shift*shift*pow(y1,1/3.0)+shift*pow(y2,1/3.0);
y[8]=shift*shift*pow(y1,1/3.0)+shift*shift*pow(y2,1/3.0);

vector<complex<double>>x;
double yx=0.0;
for(int i=0;i<9;i++)
{
yx=abs((y[i]-a[2]/3.0)*(y[i]-a[2]/3.0)*(y[i]-a[2]/3.0)+a[2]*(y[i]-a[2]/3.0)*(y[i]-a[2]/3.0)+a[1]*(y[i]-a[2]/3.0)+a[0]);
if(yx<=1.0e-13) x.push_back(y[i]-a[2]/3.0);

}
x.resize(x.size());

fstream out2("out3eq.csv", ios::out | ios::binary );
out2<<" 解（真値) "<<endl;
out2<<"x1= "<<","<<x1<<endl;
out2<<"x2= "<<","<<x2<<endl;
out2<<"x3= "<<","<<x3<<endl;
out2<<" 解（算出値) "<<endl;
for(int i=0;i<(int)x.size();i++) out2<<"x[ "<<i<<" ]= "<<x[i]<<endl;
out2.close();

remove("in3eq.csv");

return 0;
}


本3章の説明へ


1-1、三次方程式 (複素数係数) コード1-1

#include <iostream>
#include <string>
#include <cmath>
#include <complex>
#include <vector>
#include <fstream>
#include <algorithm>//transform
using namespace std;

const double pi=3.14159265358979;
const complex<double> real_unit(1.0,0.0),imag_unit(0.0,1.0);

complex<double> csum0(complex<double> d1,complex<double> d2){;return (pow(d1,1/3.0)+pow(d2,1/3.0));}

class binary_csum1
{
public: complex<double> operator()(complex<double> x,complex<double> y)
{
return pow(x,1/3.0)+pow(y,1/3.0);
}
};

int main(void)
{
complex<double> x1 = -912.0*real_unit+113.0*imag_unit, x2 = 22.0*real_unit+21.0*imag_unit,x3 = 88.0*real_unit-33.0*imag_unit;//解(真値)

ofstream out0;
out0.open("in3eq.csv", ios::out | ios::binary );
out0<<-x1*x2*x3<<endl;
out0<<(x1*x2+x1*x3+x2*x3)<<endl;
out0<<-(x1+x2+x3)<<endl;
out0.close();

vector<complex<double>>a;//a[0],a[1],a[2]となる。
complex<double> indata;

ifstream input0("in3eq.csv", ios::in | ios::binary);//
if(!input0.is_open())
{
cout<<"Cannot open file　\n";exit(1);
}
while(!input0.eof()) //eof()関数で最後までデータを読込む。
{
input0 >>indata;
if(!input0.eof()) a.push_back(indata);
}
input0.close();
a.resize(a.size());


//x^3+a[2]*x^2+a[1]*x+a[0];
complex<double> p,q;
p=a[1]-a[2]*a[2]/3.0;
q=a[0]-a[1]*a[2]/3.0+2.0*a[2]*a[2]*a[2]/27.0;


complex<double>z(0.0,0.0);
z = (q*q / 4.0 + p*p*p / 27.0);

complex<double> z1(0.0,0.0);
z1 = sqrt(z);

complex<double> y1(0.0,0.0),y2(0.0,0.0);
y1=-q/2.0+z1;
y2=-q/2.0-z1;

complex<double>c0(0.0,0.0);
transform(&y1,&y1+1,&y2,&c0,csum0);//transform使用例。

complex<double>c1(0.0,0.0);
transform(&y1,&y1+1,&y2,&c1,binary_csum1());//transform使用例。


complex<double> shift(-1.0/2.0,sqrt(3.0)/2.0);
vector<complex<double>>y(9,(0.0,0.0));//複素数解を含めると9個、必要になる。

y[0]=1.0*pow(y1,1/3.0)+1.0*pow(y2,1/3.0);
y[1]=1.0*pow(y1,1/3.0)+shift*pow(y2,1/3.0);
y[2]=1.0*pow(y1,1/3.0)+shift*shift*pow(y2,1/3.0);

y[3]=shift*pow(y1,1/3.0)+1.0*pow(y2,1/3.0);
y[4]=shift*pow(y1,1/3.0)+shift*pow(y2,1/3.0);
y[5]=shift*pow(y1,1/3.0)+shift*shift*pow(y2,1/3.0);

y[6]=shift*shift*pow(y1,1/3.0)+1.0*pow(y2,1/3.0);
y[7]=shift*shift*pow(y1,1/3.0)+shift*pow(y2,1/3.0);
y[8]=shift*shift*pow(y1,1/3.0)+shift*shift*pow(y2,1/3.0);

vector<complex<double>>x;
double yx=0.0;
for(int i=0;i<9;i++)
{
yx=abs((y[i]-a[2]/3.0)*(y[i]-a[2]/3.0)*(y[i]-a[2]/3.0)+a[2]*(y[i]-a[2]/3.0)*(y[i]-a[2]/3.0)+a[1]*(y[i]-a[2]/3.0)+a[0]);
if (yx <= 1.0e-5) x.push_back(y[i] - a[2] / 3.0);
}
x.resize(x.size());

fstream out2("out3eq.csv", ios::out | ios::binary );
out2<<" 解（真値) "<<endl;
out2<<"x1= "<<","<<x1<<endl;
out2<<"x2= "<<","<<x2<<endl;
out2<<"x3= "<<","<<x3<<endl;
out2<<" 解（算出値) "<<endl;
for(int i=0;i<(int)x.size();i++) out2<<"x[ "<<i<<" ]= ,"<<x[i]<<endl;
out2.close();

remove("in3eq.csv");

return 0;
}


本3章の説明へ



2、四次方程式 実数係数 (コード2)
#include <iostream>
#include <string>
#include <cmath>
#include <complex>
#include <vector>
#include <fstream>
#include <algorithm>
using namespace std;
const double pi=3.14159265358979;
const complex<double> real_unit(1.0,0),imag_unit(0.0,1.0);

int main(void)
{

double x10=101.1,x20=202.2,x30=303.3,x40=404.4,a40=402.0;

ofstream out0;
out0.open("in4eq.csv", ios::out | ios::binary );
out0<<x10*x20*x30*x40*a40<<endl;
out0<<-(x10*x20*x30+x10*x20*x40+x20*x30*x40+x30*x40*x10)*a40<<endl;
out0<<(x10*x20+x10*x30+x10*x40+x20*x30+x20*x40+x30*x40)*a40<<endl;
out0<<-(x10+x20+x30+x40)*a40<<endl;
out0<<a40<<endl;
out0.close();

string stringdata0;//データを文字列で読み、strtodで数値に変換し、vectorに格納する。
vector<double>a;//a[0],a[1],a[2],a[3],a[4] で、[x]内の数値は次数に対応し、係数である。
double indata;

ifstream input0("in4eq.csv", ios::in | ios::binary);
if(!input0.is_open()) {cout<<"Cannot open file　\n";exit(1);}
while(!input0.eof())
{
input0>>stringdata0;
if(!input0.eof())
{
indata=strtod((const char*)stringdata0.c_str(),NULL);//文字列を数値に変更。
a.push_back(indata);
}
}
input0.close();
a.resize(a.size());

// 3次分解式 u^3-a2/a4*u^2+b1*u+b0;
double b1,b0;
b1=a[1]*a[3]/(a[4]*a[4])-4.0*a[0]/a[4];
b0=(a[0]/a[4])*(4.0*(a[2]/a[4])-(a[3]*a[3])/(a[4]*a[4]))-(a[1]*a[1])/(a[4]*a[4]);

double p=0.0,q=0.0;
p=b1-(-a[2]/a[4])*(-a[2]/a[4])/3.0;
q=b0-b1*(-a[2]/a[4])/3.0+2.0*(-a[2]/a[4])*(-a[2]/a[4])*(-a[2]/a[4])/27.0;
complex<double>z(0.0,0.0);
z=(q*q/4.0+p*p*p/27.0)*real_unit;
complex<double> z1(0.0,0.0);
z1=pow(z,1/2.0);
complex<double> y1(0.0,0.0),y2(0.0,0.0);
y1=-q/2.0+z1;
y2=-q/2.0-z1;

complex<double> shift(-1.0/2.0,sqrt(3.0)/2.0);
vector<complex<double>>y(9,(0.0,0.0));//複素数解を含めると9個、必要になる。
y[0]=1.0*pow(y1,1.0/3.0)+1.0*pow(y2,1.0/3.0);
y[1]=1.0*pow(y1,1.0/3.0)+shift*pow(y2,1.0/3.0);
y[2]=1.0*pow(y1,1.0/3.0)+shift*shift*pow(y2,1.0/3.0);
y[3]=shift*pow(y1,1.0/3.0)+1.0*pow(y2,1.0/3.0);
y[4]=shift*pow(y1,1.0/3.0)+shift*pow(y2,1.0/3.0);
y[5]=shift*pow(y1,1.0/3.0)+shift*shift*pow(y2,1.0/3.0);
y[6]=shift*shift*pow(y1,1.0/3.0)+1.0*pow(y2,1.0/3.0);
y[7]=shift*shift*pow(y1,1.0/3.0)+shift*pow(y2,1.0/3.0);
y[8]=shift*shift*pow(y1,1.0/3.0)+shift*shift*pow(y2,1.0/3.0);

vector<double>u;
double yu=0.0;
for(int i=0;i<9;i++)
{
yu=abs((y[i]+a[2]/a[4]/3.0)*(y[i]+a[2]/a[4]/3.0)*(y[i]+a[2]/a[4]/3.0)-a[2]/a[4]*(y[i]+a[2]/a[4]/3.0)*(y[i]+a[2]/a[4]/3.0)+b1*(y[i]+a[2]/a[4]/3.0)+b0);

if(yu<=1.0e-8 && abs(imag(y[i]))<1.0e-9 && abs(imag(y[i]))/(real(y[i])+0.01)<1.0e-6) u.push_back(real(y[i])+a[2]/a[4]/3.0);
}
u.resize(u.size());

double umax=*max_element(u.begin(),u.end());

// 4次式の2次への分解 x^2+(aa-cc)x+(bb-dd) ;x^2+(aa+cc)x+(bb+dd)
double aa=a[3]/a[4]/2.0;
double bb=umax/2.0;
double dd=sqrt((umax/2.0)*(umax/2.0)-a[0]/a[4]);
double cc=0.0;

if(dd !=0.0) cc=-(a[1]/a[4]/2.0-(a[3]/a[4]/2.0)*(umax/2.0))/dd;
if(dd ==0.0) cc=sqrt(aa*aa-a[2]/a[4]+umax);

complex<double> ee1(0.0,0.0),ee2(0.0,0.0);
ee1=((aa-cc)*(aa-cc)-4.0*(bb-dd))*real_unit;
ee2=((aa+cc)*(aa+cc)-4.0*(bb+dd))*real_unit;

complex<double>x1(0.0,0.0),x2(0.0,0.0),x3(0.0,0.0),x4(0.0,0.0);

x1=-(aa-cc)/2.0*real_unit+pow(ee1,1.0/2.0)/2.0;
x2=-(aa-cc)/2.0*real_unit-pow(ee1,1.0/2.0)/2.0;
x3=-(aa+cc)/2.0*real_unit+pow(ee2,1.0/2.0)/2.0;
x4=-(aa+cc)/2.0*real_unit-pow(ee2,1.0/2.0)/2.0;

fstream out4eq("out4eq.csv", ios::out | ios::binary );
out4eq<<"解(真値)"<<endl;
out4eq<<"x10= "<<x10<<endl;
out4eq<<"x20= "<<x20<<endl;
out4eq<<"x30= "<<x30<<endl;
out4eq<<"x40= "<<x40<<endl;
out4eq<<"解(算出値) "<<endl;
out4eq<<"x1= "<<x1<<endl;
out4eq<<"x2= "<<x2<<endl;
out4eq<<"x3= "<<x3<<endl;
out4eq<<"x4= "<<x4<<endl;
out4eq.close();

remove("in4eq.csv");

return 0;
}


本3章の説明へ




2-1、四次方程式 複素数係数　コード2-1

#include <iostream>
#include <string>
#include <cmath>
#include <complex>
#include <vector>
#include <fstream>
#include <algorithm>
using namespace std;
const double pi=3.14159265358979;
const complex<double> real_unit(1.0,0),imag_unit(0.0,1.0), zero_unit(0.0, 0.0);

int main(void)
{

complex<double> x10 = 1.1*real_unit+10.0*imag_unit, x20 = 2.2*real_unit-3.0*imag_unit, x30 = -3.3*real_unit+5.0*imag_unit, x40 = 4.4*real_unit-10.0*imag_unit;


ofstream out0;
out0.open("in4eq.csv", ios::out | ios::binary );
out0<<x10*x20*x30*x40<<endl;
out0<<-(x10*x20*x30+x10*x20*x40+x20*x30*x40+x30*x40*x10)<<endl;
out0<<(x10*x20+x10*x30+x10*x40+x20*x30+x20*x40+x30*x40)<<endl;
out0<<-(x10+x20+x30+x40)<<endl;
out0<<1.0*real_unit<<endl;
out0.close();

vector<complex<double>>a;//a[0],a[1],a[2],a[3],a[4] で、[x]内の数値は次数に対応し、係数である。
complex<double> indata;

ifstream input0("in4eq.csv", ios::in | ios::binary);
if(!input0.is_open()) {cout<<"Cannot open file　\n";exit(1);}
while(!input0.eof())
{
input0 >> indata;
if(!input0.eof()) a.push_back(indata);
}
input0.close();
a.resize(a.size());


// 3次分解式 u^3-a2/a4*u^2+b1*u+b0;
complex<double> b1, b0;

b1=a[1]*a[3]/(a[4]*a[4])-4.0*a[0]/a[4];
b0=(a[0]/a[4])*(4.0*(a[2]/a[4])-(a[3]*a[3])/(a[4]*a[4]))-(a[1]*a[1])/(a[4]*a[4]);


complex<double> p = zero_unit, q = zero_unit;

p=b1-(-a[2]/a[4])*(-a[2]/a[4])/3.0;
q=b0-b1*(-a[2]/a[4])/3.0+2.0*(-a[2]/a[4])*(-a[2]/a[4])*(-a[2]/a[4])/27.0;

complex<double>z(0.0, 0.0);
z = (q*q / 4.0 + p*p*p / 27.0);

complex<double> z1(0.0,0.0);
z1=pow(z,1/2.0);

complex<double> y1(0.0,0.0),y2(0.0,0.0);
y1=-q/2.0+z1;
y2=-q/2.0-z1;

complex<double> shift(-1.0/2.0,sqrt(3.0)/2.0);
vector<complex<double>>y(9,(0.0,0.0));//複素数解を含めると9個、必要になる。
y[0]=1.0*pow(y1,1.0/3.0)+1.0*pow(y2,1.0/3.0);
y[1]=1.0*pow(y1,1.0/3.0)+shift*pow(y2,1.0/3.0);
y[2]=1.0*pow(y1,1.0/3.0)+shift*shift*pow(y2,1.0/3.0);
y[3]=shift*pow(y1,1.0/3.0)+1.0*pow(y2,1.0/3.0);
y[4]=shift*pow(y1,1.0/3.0)+shift*pow(y2,1.0/3.0);
y[5]=shift*pow(y1,1.0/3.0)+shift*shift*pow(y2,1.0/3.0);
y[6]=shift*shift*pow(y1,1.0/3.0)+1.0*pow(y2,1.0/3.0);
y[7]=shift*shift*pow(y1,1.0/3.0)+shift*pow(y2,1.0/3.0);
y[8]=shift*shift*pow(y1,1.0/3.0)+shift*shift*pow(y2,1.0/3.0);


vector<complex<double>>u;
double yu=0.0;
for(int i=0;i<9;i++)
{
yu=abs((y[i]+a[2]/a[4]/3.0)*(y[i]+a[2]/a[4]/3.0)*(y[i]+a[2]/a[4]/3.0)-a[2]/a[4]*(y[i]+a[2]/a[4]/3.0)*(y[i]+a[2]/a[4]/3.0)+b1*(y[i]+a[2]/a[4]/3.0)+b0);

if (yu <= 1.0e-5 ) u.push_back(y[i] + a[2] / a[4] / 3.0);
}
u.resize(u.size());

vector<double> uabs(u.size(),0.0);
for (int i = 0; i <(int)u.size(); i++) uabs[i] = abs(u[i]);
int uabs_max=max_element(uabs.begin(),uabs.end())-uabs.begin();//絶対最大値の位置探しで uabs.begin()で引き算している。
complex<double> umax =u[uabs_max];//絶対値最大の複素数解


// 4次式の2次への分解 x^2+(aa-cc)x+(bb-dd) ;x^2+(aa+cc)x+(bb+dd)
complex<double> aa = a[3] / a[4] / 2.0;
complex<double> bb = umax / 2.0;
complex<double> dd = sqrt((umax / 2.0)*(umax / 2.0) - a[0] / a[4]);
complex<double> cc = zero_unit;

if (dd != zero_unit) cc = -(a[1] / a[4] / 2.0 - (a[3] / a[4] / 2.0)*(umax / 2.0)) / dd;
if (dd == zero_unit) cc = sqrt(aa*aa - a[2] / a[4] + umax);

complex<double> ee1(0.0,0.0),ee2(0.0,0.0);
ee1 = ((aa - cc)*(aa - cc) - 4.0*(bb - dd));
ee2 = ((aa + cc)*(aa + cc) - 4.0*(bb + dd));

complex<double>x1(0.0,0.0),x2(0.0,0.0),x3(0.0,0.0),x4(0.0,0.0);

x1 = -(aa - cc) / 2.0 + pow(ee1, 1.0 / 2.0) / 2.0;
x2 = -(aa - cc) / 2.0 - pow(ee1, 1.0 / 2.0) / 2.0;
x3 = -(aa + cc) / 2.0 + pow(ee2, 1.0 / 2.0) / 2.0;
x4 = -(aa + cc) / 2.0 - pow(ee2, 1.0 / 2.0) / 2.0;

fstream out4eq("out4eq.csv", ios::out | ios::binary );
out4eq<<"解(真値)"<<endl;
out4eq<<"x10=, "<<x10<<endl;
out4eq<<"x20=, "<<x20<<endl;
out4eq<<"x30=, "<<x30<<endl;
out4eq<<"x40=, "<<x40<<endl;
out4eq<<"解(算出値) "<<endl;
out4eq<<"x1=, "<<x1<<endl;
out4eq<<"x2=, "<<x2<<endl;
out4eq<<"x3=, "<<x3<<endl;
out4eq<<"x4=, "<<x4<<endl;
out4eq.close();

remove("in4eq.csv");

return 0;
}


本3章の説明へ



3、Newton-Raphson法による方程式の求根(コード3)
#include <cmath>
#include <fstream>
using namespace std;

class f_orig { public: double operator()(double x) { return x*x*x-5.2*x*x+2.2*x+32.2;} };
//解 -1.9711631562230,他2根は複素数

class f_diff { public: double operator()(double x) { return 3.0*x*x-5.2*2.0*x+2.2;} };

int main() 　
{
double eps=1.0e-10;
double a0=1.0;
double a1=-5.2;
double x=-a1/a0;

while(fabs(f_orig()(x))>eps) x=x-f_orig()(x)/f_diff()(x);

ofstream cout0("outnewton.txt",ios::out |ios::binary);
cout0<<" x= "<<x<<endl;
cout0.close();

return 0;
}


本3章の説明へ


4、二分法類似法(コード4)
#include <cmath>
#include <fstream>
#include <vector>
#include <algorithm>//find_if、transform
#include <functional>// multiplies
using namespace std;

class f_bisection { public: double operator()(double x) { return x/2.0*exp(x*x)-exp(x); } };

bool minusvalue(double xx) { if(xx<0.0) return true; return false; }

//区間細分法（二分法の一種と考える。）
//区間を細分した時、ある区間内にゼロ点あるとその前後で、f(x)の値が変化する。
//区分が一つ増加部との積、v[n]*v[n+1]は必ず負値になる。その区間以外は全て正である。
//この考え方で、zero点区間をSTLを利用し求めて、解を算出する。

int main()
{
double eps=1.0e-12;
double a,b;

b=5.0;//解を知らない時は、カットアンドトライで望むしかない。
a=1.0;

int n=(int)((log((b-a)/eps)/log(2.0)+0.5)/3.00);
double h=(b-a)/n;

vector<double>v1(n+1,0.0);//　n+1の数は位置bでのf(b)値を含めるため。
for(int i=0;i<=n;i++) v1[i]=f_bisection()(a+i*h);

vector<double>v2(n+1,0.0); //f(a+(n-1)*h)*f(b)の積を求めるため。
transform(v1.begin(),v1.end()-1,v1.begin()+1,v2.begin(),multiplies<double>());//積をv2に代入。


// 負値になる区間を算出する。
vector<double>::iterator p1;// find_if()のitaratorである。
p1=find_if(v2.begin(),v2.end(),minusvalue);//p1で負値の個所を得る。
int near0=p1-v2.begin();


int j=0;//収束のための繰返し計算
double solv;

for(int k=0;k<20;k++)
{
b=a+(near0+1)*h;
a=a+near0*h;
h=(b-a)/n;
j++;
vector<double>v3(n+1,0.0);
for(int i=0;i<=n;i++) v3[i]=f_bisection()(a+i*h);
vector<double>v4(n+1,0.0);
transform(v3.begin(),v3.end()-1,v3.begin()+1,v4.begin(),multiplies<double>());
int near01=find_if(v4.begin(),v4.end(),minusvalue)-v4.begin();//簡略文

near0=near01;
v3.clear();
v4.clear();

solv=(a+near01*h+a+(near01+1)*h)/2;
if(fabs(f_bisection()(solv))<eps) break;
}

ofstream cout35("calculation_results.csv",ios::out |ios::binary );
cout35<<" calculation results "<<endl;
cout35.precision(15);
cout35<<" solved value = "<<solv<<endl;
cout35<<" "<<endl;
cout35<<" 計算繰り返し数 = "<<j<<endl;
cout35<<" x=solved value、f(x) = "<<f_bisection()(solv)<<endl;
cout35.close();

return 0;
}


本3章の説明へ


5、常微分方程式 Runge-Kutta(コード5)
#include <cmath>
#include <fstream>
#include <vector>
#include <functional>//binary_function
using namespace std;
const double pi=3.1415926535897932;

class function_xy { public: double operator()(double x, double y) { return x+x*cos(y);} };

struct f_xy:public binary_function<double,double,double>{ double operator()(double x,double y) {return x+x*cos(y);} };

// 4次のrunnge-kutta法による微分方程式の解法
// typedif文とclass文のfunction文としての使用法を示した。
int main()
{

double x,y;
double k,k1,k2,k3,k4;
int n=60;
double h=pi/(float)n;
int i=0;
x=y=0;

//class文の使用例
ofstream cout1("out_class.csv",ios::out |ios::binary);
cout1<<" i "<<","<<" x "<<","<<" y "<<endl;
while(i<n)
{
i++;
k1=h*function_xy()(x,y);//class文の使用例
k2=h*function_xy()(x+h/2.0,y+k1/2.0);
k3=h*function_xy()(x+h/2.0,y+k2/2.0);
k4=h*function_xy()(x+h,y+k3);
k=(k1+2.0*k2+2.0*k3+k4)/6.0;
y=y+k;
x=x+h;
cout1<<i<<","<<x<<","<<y<<endl;
}
cout1.close();

//typedif文の使用例
vector<double>fx;
vector<double>fy;
i=0;
x=y=0;
while(i<n)
{
i++;
k1=h*f_xy()(x,y);
k2=h*f_xy()(x+h/2.0,y+k1/2.0);
k3=h*f_xy()(x+h/2.0,y+k2/2.0);
k4=h*f_xy()(x+h,y+k3);
k=(k1+2.0*k2+2.0*k3+k4)/6.0;
y=y+k;
x=x+h;
fy.push_back(y);
fx.push_back(x);
}
fy.resize(fy.size());
fx.resize(fx.size());

ofstream cout2("out_struct.csv",ios::out |ios::binary);
cout2<<" i "<<","<<" x "<<","<<" y "<<endl;
for(int i=0;i<(int)fy.size();++i) cout2<<i<<","<<fx[i]<<","<<fy[i]<<endl;
cout2.close();

return 0;
}


本3章の説明へ


6、最小二乗法　直線近似(コード6)
#include <cmath>
#include <fstream>
#include <vector>
#include <numeric> // for inner_product
using namespace std;

int main()
{

//1、最小二乗法による直線近似 近似のa*x+bのaとbの算出

int nn=10; //データ数
vector<double>u;
u.assign(2*nn,0.0);

// y=a*x+b=3.6*x+0.3
// x軸
u[0]=3.0;
u[1]=3.3;
u[2]=3.6;
u[3]=4.0;
u[4]=4.1;
u[5]=5.2;
u[6]=5.4;
u[7]=6.0;
u[8]=7.0;
u[9]=8.6;

// y軸 下記の付加している数値(例 +1.5)は最小二乗の計算のために、正常な値に追加した数値である。
u[10]=11.1;
u[11]=12.18+1.5;
u[12]=13.26;
u[13]=14.7;
u[14]=15.06+2.0;
u[15]=19.02;
u[16]=19.74-3.0;
u[17]=21.9;
u[18]=25.5-4.0;
u[19]=31.26+5.0;

//u.size()/2は、データ数nnと同じ。

double average_x=accumulate(u.begin(),u.begin()+u.size()/2,0.0)/(u.size()/2);// xの平均
double average_y=accumulate(u.begin()+u.size()/2,u.end(),0.0)/(u.size()/2);// yの平均
double average_xy=inner_product(u.begin(),u.begin()+u.size()/2,u.begin()+u.size()/2,0.0)/(u.size()/2);//ｘとｙの積の平均　　　
double average_xx=inner_product(u.begin(),u.begin()+u.size()/2,u.begin(),0.0)/(u.size()/2);//xの二乗の平均
double a=(average_xy-average_x*average_y)/(average_xx-average_x*average_x);
double b=(average_y*average_xx-average_x*average_xy)/(average_xx-average_x*average_x);

ofstream coutw("least_linear.csv",ios::out |ios::binary);
coutw<<"a"<<","<<"b"<<endl;
coutw<<a<<","<<b<<endl;
coutw.close();

//2、相関係数R^2の算出（残差、偏差の2法で算出）
//相関係数 the Pearson product-moment correlation coefficient (the PMCC）でrで表現。

// 1)ε残差とｙ残差による相関係数算出
//a)STL使用　相関係数をrsq1a_pmcc。 　　
double sum_res2a=0.0;
for(int i=0;i<(int)u.size()/2;i++) sum_res2a+=(u[i+u.size()/2]-(a*u[i]+b))*(u[i+u.size()/2]-(a*u[i]+b));
double syy=0.0;
for(int i=0;i<(int)u.size()/2;i++) syy+=(u[i+u.size()/2]-average_y)*(u[i+u.size()/2]-average_y);

double rsq1a_pmcc=1-sum_res2a/syy;

//b）STLの極力使用 相関係数をrsq1b_pmcc
vector<double>vaub;
for(int i=0;i<(int)u.size()/2;i++) vaub.push_back(u[i+u.size()/2]-(a*u[i]+b));
vaub.resize(vaub.size());
double sum_res2b=inner_product(vaub.begin(),vaub.end(),vaub.begin(),0.0);
vector<double>vyy;
for(int i=0;i<(int)u.size()/2;i++) vyy.push_back(u[i+u.size()/2]-average_y);
vyy.resize(vyy.size());
double tyy=inner_product(vyy.begin(),vyy.end(),vyy.begin(),0.0);

double rsq1b_pmcc=1-sum_res2b/tyy;

// 2)xyの偏差値による算出 相関係数をrsqc_pmcc
double s_resx1=0.0;
double s_resy1=0.0;
double s_resxy1=0.0;
for(int i=0;i<(int)u.size()/2;i++) s_resx1+=(u[i]-average_x)*(u[i]-average_x);
for(int i=0;i<(int)u.size()/2;i++) s_resy1+=(u[i+u.size()/2]-average_y)*(u[i+u.size()/2]-average_y);
for(int i=0;i<(int)u.size()/2;i++) s_resxy1+=(u[i]-average_x)*(u[i+u.size()/2]-average_y);

double rsqc_pmcc=(s_resxy1)*(s_resxy1)/s_resx1/s_resy1;

ofstream coutw1("pmcc_linear.csv",ios::out |ios::binary|ios::app);
coutw1<<"rsq1a_pmcc= "<<","<<rsq1a_pmcc<<endl;
coutw1<<"rsq1b_pmcc= "<<","<<rsq1b_pmcc<<endl;
coutw1<<"rsqc_pmcc= "<<","<<rsqc_pmcc<<endl;
coutw1.close();

return 0;
}


本3章の説明へ



7、最小二乗 多項式近似 データ列数1列(コード7)
#include <string>
#include <cmath>
#include <fstream>
#include <vector>
#include <algorithm>
#include <iomanip>
#include <numeric>
using namespace std;

class nthpower { public: double operator()(double sx, double sn) { return pow(sx,sn); } };

int main()
{
// 最小二乗法　多項式近似 多項式近似値算出には、Gauss-Jordannを使用する。 相関係数(pmcc)のr-2を残余から求める。

//0、近似用データの生成
int mm=20;//データ数
vector<double>xx;// x軸のデータ用
vector<double>yy;//　y=f(x)のy軸のデータ用
xx.assign(mm,0.0);//初期化
yy.assign(mm,0.0);

for(int i=0;i<mm;++i) xx[i]=0.50*i+i*i/500.0+5.00;
for(int i=0;i<mm;++i) yy[i]=+50.0+5.0*xx[i]-10.0*xx[i]*xx[i]+0.5*xx[i]*xx[i]*xx[i]+0.05*xx[i]*xx[i]*xx[i]*xx[i];
yy[1]=yy[1]+15.6;
yy[3]=yy[3]+22.8;
yy[6]=yy[6]-13.6;
yy[8]=yy[8]+20.6;

ofstream outxxyy("sample_data.csv",ios::out |ios::binary );//データfile
outxxyy<<"データ名"<<endl;
outxxyy<<"近似次数"<<endl;
outxxyy<<4<<endl;
outxxyy<<"データ番号"<<","<<"x値"<<","<<"y値"<<endl;
for(int i=0;i<mm;i++) outxxyy<<i<<","<<xx[i]<<","<< yy[i]<<endl;
outxxyy.close();


//1、データfileからデータの読み込み、// データは上記のようにデータ番号、x、yと横3列とする。
//1-1、 STL/stringで読み込み

string string0,string1,string2;
string stringdata0;
vector<string>stringdata;
int npmax;

ifstream input0("sample_data.csv", ios::in | ios::binary);
if(!input0.is_open()) exit(1);
input0>>string0;
input0>>string1;
input0>>npmax;
input0>>string2;
while(!input0.eof())
{
input0>>stringdata0;
if(!input0.eof())
{
stringdata.push_back(stringdata0);
}
}
input0.close();
stringdata.resize(stringdata.size());


//1-2、string型dataである読込んだstringdata1の数値データ化 //stringdata[j] 列の数は3に固定
vector<double>xv;//数値入力データの格納
string conma=",";// note;下記 "place" is "size_type";

for(int j=0;j<(int)stringdata.size();j++)
{
string ::iterator place0=find_first_of(stringdata[j].begin(),stringdata[j].end(),conma.begin(),conma.end());//<algorithm>
if(place0 ==stringdata[j].end()) break;

string s0;//文字列数値の切り出し
for(int i=0;i<(place0-stringdata[j].begin());i++) s0.insert(s0.end(),stringdata[j][i]); //stringdata[j][i]ができる。

//文字列の数値化
string s0charcheck=s0;
const char* check0=s0charcheck.c_str();
double xs0;
xs0 = strtod(check0, NULL);
xv.push_back(xs0);

// 取りだした文字列を削除していく。tecnique上”+1”が必要。
basic_string <char>::iterator strg0_Iter;
strg0_Iter =stringdata[j].erase (stringdata[j].begin(),stringdata[j].begin()+(place0-stringdata[j].begin())+1 );

//新規のvectorとなったstringdata[j]の文字列から","までの位置を求める。
string::iterator place1=find_first_of(stringdata[j].begin(),stringdata[j].end(),conma.begin(),conma.end());
string s1;
for(int i=0;i<(place1-stringdata[j].begin());i++) s1.insert(s1.end(),stringdata[j][i]);
string s1charcheck=s1;
const char* check1=s1charcheck.c_str();
double xs1;
xs1 = strtod(check1,NULL);
xv.push_back(xs1);
basic_string <char>::iterator strg1_Iter;
strg1_Iter =stringdata[j].erase (stringdata[j].begin(),stringdata[j].begin()+(place1-stringdata[j].begin())+1 );

string::iterator place2=find_first_of(stringdata[j].begin(),stringdata[j].end(),conma.begin(),conma.end());
string s2;
for(int i=0;i<(place2-stringdata[j].begin());i++) s2.insert(s2.end(),stringdata[j][i]);
string s2charcheck=s2;
const char* check2=s2charcheck.c_str();
double xs2;
xs2 = strtod(check2, NULL);
xv.push_back(xs2);
basic_string <char>::iterator strg2_Iter;
strg2_Iter=stringdata[j].erase(stringdata[j].begin(),stringdata[j].begin()+(place2-stringdata[j].begin()));//+1が無い事に注意。

}
xv.resize(xv.size());

vector<double>x,y;// 3で割っているのは列が3であるため
for(int i=0;i<(int)xv.size()/3;i++) //xvの最初がデータ番号である事が解っているので、
{
x.push_back(xv[3*i+1]);// 3*iの3は3個毎に同様データが繰り返されるので
y.push_back(xv[3*i+2]);
}
x.resize(x.size());
y.resize(y.size());


// 2、最小二乗法　多項式近似 次数np
// 多項式近似値算出には、Gauss-Jordannの連立一次方程式の解を算出を使用する。相関係数(pmcc)のr-2を残余から求める方法で算出する。

//2-1、データ数をm、求める多項式の次数をnp Σxiの(i<m)の一次、二次、三次、…、2(np+1)次までの係数を作成する。(np+1)(np+1)正方行列のデータを解く。


int np;//npは次数。データfileからの入力npmaxで決める。
np=npmax;
int m=(int)xv.size()/3;

vector<double>v1,v2;
v1.assign(m,0.0);
v2.push_back(m);//多項式近似の計算よりデータ数mを最初に。push_frontはNO。v2.assign(2m,0)ではないのは、push_back値が2mの後となる。

//2-2、Σxiの2*np乗のデータ作成
vector<int>norder;

//1）x側の値
double sum_mdata;
for(int i=1;i<=2*np;++i)
{
norder.assign(m,i);
transform(x.begin(),x.end(),norder.begin(),v1.begin(),nthpower());//xiのi乗をnthpowerで生成後、v1に代入。
sum_mdata=accumulate(v1.begin(),v1.end(),0.0);
v2.push_back(sum_mdata);
}

// 2）y側の値
vector<double>v3,v4;
v3.assign(m,0.0);
v4.push_back(accumulate(y.begin(),y.end(),0.0));
double resultxy;
for(int i=1;i<=2*np;++i)
{
norder.assign(m,i);
transform(x.begin(),x.end(),norder.begin(),v3.begin(),nthpower());
resultxy=inner_product(v3.begin(),v3.end(),y.begin(),0.0);
v4.push_back(resultxy);
}


//2-3、連立方程式の適用
vector<double>s,t;//連立方程式の係数側 　
vector<double>ans;// 解を収納するベクトル
vector<double>ddd;
vector<int>row;
s.clear();
t.clear();
ans.clear();
ddd.clear();
row.clear();

//2-4、連立方程式の係数作成
int i,j,k;
int n=np+1;
t.assign(n,0.0);
for(int j=0;j<n;j++)
{
for(int i=j;i<n+j;i++) s.push_back(v2[i]);
}
s.resize(s.size());
for(int i=0;i<n;i++) t[i]=v4[i];

vector<double>ss,tt;
ss.assign(n*n,0.0);
tt.assign(n,0.0);
ss=s;
tt=t;
ans.assign(n,0.0);
for(i=0;i<n;++i) row.insert(row.end(),1,i);
for(k=0;k<n;k++)
{
ddd.assign(n*n,0.0);
transform(s.begin()+n*k,s.end(),ddd.begin(),fabs);

int inum=max_element(ddd.begin(),ddd.end())-ddd.begin();
int maxvline=(int)((n*k+inum)/n);
int maxvrow=n*k+inum- maxvline*n;
for(i=0;i<n;i++) swap(s[n*k+i],s[n*(maxvline)+i]);
for(i=0;i<n;i++) swap(s[i*n+k],s[i*n+maxvrow]);
swap(t[k],t[maxvline]);
swap(row[k],row[maxvrow]);
for (i=k+1; i<n; ++i)
{
t[i]-=t[k]*s[i*n+k]/s[k*n+k];
for(j=n-1;j>=k;j--)
{
if(fabs(s[k*n+k])>1.0e-15) s[i*n+j] -=s[i*n+k]*s[k*n+j]/s[k*n+k];
}
}
}
double anssum;
anssum=0.0;
for(i=n-1;i>=0;--i)
{
anssum=0.0;
for(j=0;j<n-i-1;j++) anssum=s[n*i+n-1-j]*ans[row[n-1-j]]+anssum;
ans[row[i]]=(t[i]-anssum)/s[n*i+i];
}


//2-5、計算結果検証用データ作成
vector<double>recal;
double sum;
for(j=0;j<m;j++)
{
sum=0.0;
for(i=1;i<n;i++) sum=sum+ans[i]*pow(x[j],i);
recal.push_back(sum+ans[0]);
}
recal.resize(recal.size());

//2-5、相関係数R^2の算出
// 1)ε残差とｙ残差による相関係数算出
double average_x=accumulate(x.begin(),x.end(),0.0)/(x.size());
double average_y=accumulate(y.begin(),y.end(),0.0)/(y.size());
double sum_resty=0.0;
for(int i=0;i<(int)y.size();i++) sum_resty+=(y[i]-recal[i])*(y[i]-recal[i]);
double sum_yy=0.0;
for(int i=0;i<(int)y.size();i++) sum_yy+=(y[i]-average_y)*(y[i]-average_y);
double rsqrest_pmcc=1-sum_resty/sum_yy;

ofstream out1("Calculated_value.csv",ios::out |ios::binary);
out1<<string0<<endl;
out1<<string1<<endl;
out1<<"多項式 入力次数= "<<npmax<<endl;
out1<<" "<<endl;
out1<<"計算結果 近似係数の順＝定数、1次…n次　"<<endl;
for(i=0;i<(int)ans.size();i++)
out1<<setw(2)<<i<<","<<setw(15)<<setprecision(10)<<ans[i]<<endl;
out1<<"近似の検証　"<<endl;
out1<<" i "<<","<<" re-cal_value "<<","<<" y[i] "<<","<<"recal-y[i]"<<endl;
for(i=0;i<(int)recal.size();i++)
out1<<i<<","<<recal[i]<<","<< y[i]<<","<<recal[i]-y[i]<<endl;
out1<<" "<<endl;
out1<<"相関係数R2 "<<","<<rsqrest_pmcc<<endl;
out1.close();

remove("sample_data.csv");

return 0;
}


本3章の説明へ



8、最小二乗 多項式近似 データ列自由(コード8)
#include <string>
#include <cmath>
#include <fstream>
#include <vector>
#include <algorithm>
#include <functional>
#include <iomanip>
#include <numeric>
using namespace std;

class nthmultiplies
{
public:
double operator()( double x, double y)
{
int ny=(int)y;
double mul0=1.0;
for(int i=0;i<ny;i++) mul0*=x;
return mul0;
}
};


int main()
{
// 最小二乗法　多項式近似 Gauss-Jordanの解を算出を使用。 相関係数(pmcc)のr-2を残余から求める方法で算出。

//0、近似用データの生成
int mm=50;//データ数
vector<double>xx;//xのデータ格納
vector<double>yy0;//y0=f(x)のy0のデータ
vector<double>yy1;
vector<double>yy2;
vector<double>yy3;
vector<double>yy4;
xx.assign(mm,0.0);
yy0.assign(mm,0.0);
yy1.assign(mm,0.0);
yy2.assign(mm,0.0);
yy3.assign(mm,0.0);
yy4.assign(mm,0.0);
double xt=0.0;

for(int i=0;i<mm;++i)
{
xt=((int)(10.0*(0.50*i+i*i/500.0-10.00)))/10.0;
xx[i]=xt;
yy0[i]=5.0+0.50*xt-2.0*xt*xt;
yy1[i]=5.0+0.5*xt-2.0*xt*xt+1.0*xt*xt*xt;
yy2[i]=5.0+0.5*xt-2.0*xt*xt+1.0*xt*xt*xt-1.0*xt*xt*xt*xt;
yy3[i]=5.0+0.5*xt-2.0*xt*xt+1.0*xt*xt*xt-1.0*xt*xt*xt*xt+0.5*xt*xt*xt*xt*xt;
yy4[i]=5.0+0.5*xt-2.0*xt*xt+1.0*xt*xt*xt-1.0*xt*xt*xt*xt+0.5*xt*xt*xt*xt*xt+1.0*xt*xt*xt*xt*xt*xt;
}
ofstream outxxyy("sample_data.csv",ios::out |ios::binary );//データfile
outxxyy<<"データ名sample兼入力データ作成方法 "<<endl;
outxxyy<<"近似次数最大値"<<endl;
outxxyy<<8<<endl;
outxxyy<<"x値とする列番号の指定(仮設定1…0起点故2列目で1)"<<endl;
outxxyy<<1<<endl;
outxxyy<<"近似対象データの列並び(仮設定2…0起点故3列目以右で2)"<<endl;
outxxyy<<2<<endl;
outxxyy<<"データ番号"<<","<<"x値"<<","<<"y値(以右列)"<<endl;
for(int i=0;i<mm;i++)
outxxyy<<i<<","<<xx[i]<<","<< yy0[i]<<","<<yy1[i]<<","<<yy2[i]<<","<<yy3[i]<<","<<yy4[i]<<endl;
outxxyy.close();

//1、csvデータの読み込み（STL/string使用）
// データ列に制限はない。（上記例は、x1列、y横5列であるが、計算可能）
// データ行数に制限はない。

string string0,string1,string2,string3,string4;
string stringdata0;
vector<string>stringdata;
int npmax;
int x_location;
int y_location;
ifstream input0("sample_data.csv", ios::in | ios::binary);
if(!input0.is_open()) exit(1);
input0>>string0;
input0>>string1;
input0>>npmax;
input0>>string2;
input0>>x_location;
input0>>string3;
input0>>y_location;
input0>>string4;

while(!input0.eof())
{
input0>>stringdata0;
if(!input0.eof())
{
stringdata.push_back(stringdata0);
}
}
input0.close();
stringdata.resize(stringdata.size());

//2、string型dataであるstringdata1の数値データ化 列数算出　stringdata[j]の末尾に","を付け加える。
for(int j=0;j<(int)stringdata.size();j++) stringdata[j]=stringdata[j]+",";

vector<double>xv;
string conma=",";
int kp;//csv列数
for(int j=0;j<(int)stringdata.size();j++)
{
kp=0;
while(stringdata[j].size())
{
string ::iterator place0=find_first_of(stringdata[j].begin(),stringdata[j].end(),conma.begin(),conma.end());
++kp;
if(place0 !=stringdata[j].end())
{
string s0;
for(int i=0;i<(place0-stringdata[j].begin());i++) s0.insert(s0.end(),stringdata[j][i]); //stringdata[j][i]可能。
string s0charcheck=s0;
const char* check0=s0charcheck.c_str();
double xs0;
xs0 = strtod(check0, NULL);//文字列の数値化
xv.push_back(xs0);
basic_string <char>::iterator strg0_Iter;
strg0_Iter =stringdata[j].erase (stringdata[j].begin(),stringdata[j].begin()+(place0-stringdata[j].begin())+1);
}
else if(place0==(stringdata[j].end()))
{
string s1;
for(int i=0;i<(place0-stringdata[j].begin());i++) s1.insert(s1.end(),stringdata[j][i]);
string s1charcheck=s1;
const char* check1=s1charcheck.c_str();
double xs1;
xs1 = strtod(check1, NULL);
xv.push_back(xs1);
basic_string <char>::iterator strg0_Iter;
strg0_Iter =stringdata[j].erase (stringdata[j].begin(),stringdata[j].begin()+(place0-stringdata[j].begin()));
break;
}
}
}
xv.resize(xv.size());

//3、xとｙデータ列の組み合わせ // x列は[kp*i+x_location]で変更可 // y列はnum_rowでy_locationから最後まで計算
vector<double>x,y;
for(int num_row=y_location;num_row<kp;num_row++)//return直前の" ｝"と対応
{

x.clear();
y.clear();

for(int i=0;i<(int)xv.size()/kp;i++)
{
x.push_back(xv[kp*i+x_location]);
y.push_back(xv[kp*i+num_row]);//ｙ列の選定
}
x.resize(x.size());
y.resize(y.size());

//4、最小二乗法　多項式近似値算出:Gauss-Jordan 相関係数(pmcc)のr-2:残余から得る。
// データ数をm、求める多項式の次数をnp
// Σxiの(i<m)の一次、二次、三次、…、2(np+1)次までの係数を作成する。
// 最終的には、(np+1)(np+1)正方行列のデータになる。

//4-1、Σxiの2*np乗のデータ作成
int np;
int m=(int)xv.size()/kp;
vector<double>v1,v2;
vector<int>norder;
vector<double>v3,v4;
for(np=2;np<=npmax;np++) //近似の次数を2からnpmaxを実施する。return二つ前の" ｝"と対応
{
norder.clear();
v1.clear();
v2.clear();
v1.assign(m,0.0);
v2.push_back(m);//データ数mを最初の位置に入れる。

//4-2、2*np乗値算出用transfomer(2項関数(nthmultiplies()))による正方行列用データ作成
//1）x側の値
double sum_mdata=0.0;
for(int i=1;i<=2*np;++i)
{
norder.assign(m,i);
transform(x.begin(),x.end(),norder.begin(),v1.begin(),nthmultiplies());
sum_mdata=accumulate(v1.begin(),v1.end(),0.0);
v2.push_back(sum_mdata);
}
// 2）y側の値
v4.clear();
v3.assign(m,0.0);
v4.push_back(accumulate(y.begin(),y.end(),0.0));
double resultxy=0.0;
for(int i=1;i<=2*np;++i)
{
norder.assign(m,i);
transform(x.begin(),x.end(),norder.begin(),v3.begin(),nthmultiplies());
resultxy=inner_product(v3.begin(),v3.end(),y.begin(),0.0);
v4.push_back(resultxy);
}

//4-3、連立方程式
vector<double>s,t;//連立方程式の係数側 　
vector<double>ans;//解ベクトル
vector<double>ddd;
vector<int>row;
s.clear();
t.clear();
ans.clear();
ddd.clear();
row.clear();

//4-4、連立方程式用の係数作成（連立方程式は、ｎ次項に対して(np+1)(np+1)の正方行列）
int i,j,k;
int n=np+1;
t.assign(n,0.0);
for(int j=0;j<n;j++)
{
for(int i=j;i<n+j;i++) s.push_back(v2[i]);
}
s.resize(s.size());
for(int i=0;i<n;i++) t[i]=v4[i];

vector<double>ss,tt;
ss.assign(n*n,0.0);
tt.assign(n,0.0);
ss=s;
tt=t;
ans.assign(n,0.0);
for(i=0;i<n;++i) row.insert(row.end(),1,i);
for(k=0;k<n;k++)
{
ddd.assign(n*n,0.0);
transform(s.begin()+n*k,s.end(),ddd.begin(),fabs);
int inum=max_element(ddd.begin(),ddd.end())-ddd.begin();
int maxvline=(int)((n*k+inum)/n);
int maxvrow=n*k+inum- maxvline*n;
for(i=0;i<n;i++) swap(s[n*k+i],s[n*(maxvline)+i]);
for(i=0;i<n;i++) swap(s[i*n+k],s[i*n+maxvrow]);
swap(t[k],t[maxvline]);
swap(row[k],row[maxvrow]);
for (i=k+1; i<n; ++i)
{
t[i]-=t[k]*s[i*n+k]/s[k*n+k];
for(j=n-1;j>=k;j--)
{
if(fabs(s[k*n+k])>1.0e-15) s[i*n+j] -=s[i*n+k]*s[k*n+j]/s[k*n+k];
}
}
}
double anssum;
anssum=0.0;
for(i=n-1;i>=0;--i)
{
anssum=0.0;
for(j=0;j<n-i-1;j++) anssum=s[n*i+n-1-j]*ans[row[n-1-j]]+anssum;
ans[row[i]]=(t[i]-anssum)/s[n*i+i];
}

//5、計算結果検証
vector<double>recal;
double sum;
for(j=0;j<m;j++)
{
sum=0.0;
for(i=1;i<n;i++) sum=sum+ans[i]*pow(x[j],i);
recal.push_back(sum+ans[0]);
}
recal.resize(recal.size());

//6、相関係数R^2の算出
// 1)ε残差とｙ残差による相関係数算出
double average_x=accumulate(x.begin(),x.end(),0.0)/(x.size());
double average_y=accumulate(y.begin(),y.end(),0.0)/(y.size());
double sum_resty=0.0;
for(int i=0;i<(int)y.size();i++) sum_resty+=(y[i]-recal[i])*(y[i]-recal[i]);
double sum_yy=0.0;
for(int i=0;i<(int)y.size();i++) sum_yy+=(y[i]-average_y)*(y[i]-average_y);
double rsqrest_pmcc=1-sum_resty/sum_yy;

ofstream out1("Calculation_result.csv",ios::out |ios::binary |ios::app);
out1<<" "<<endl;
out1<<string0<<endl;
out1<<"指定した計算の最高次数(入力値)= "<<npmax<<endl;
out1<<"多項式　x値の列番号= "<<x_location<<" 近似対象データ列番号= "<<num_row<<endl;
out1<<"多項式適用次数= "<<np<<endl;
out1<<" "<<endl;
out1<<"計算結果 近似係数の順＝定数、1次…n次　"<<endl;
out1<<"次数"<<","<<"近似項の係数"<<endl;
for(i=0;i<(int)ans.size();i++)
out1<<setw(2)<<i<<","<<setw(15)<<setprecision(10)<<ans[i]<<endl;
out1<<" "<<endl;
out1<<"近似式の検証　"<<endl;
out1<<" i "<<","<<"x値"<<","<<" 再計算値（近似式）"<<","<<" 原値(y)"<<","<<"再計算値-原値"<<endl;
for(i=0;i<(int)recal.size();i++)
out1<<i<<","<<x[i]<<","<<recal[i]<<","<< y[i]<<","<<recal[i]-y[i]<<endl;
out1<<" "<<endl;
out1<<"相関係数R2 "<<","<<rsqrest_pmcc<<endl;
out1<<" "<<endl;
out1<<" "<<endl;
out1.close();

}// 多項式次数npのforに対応
}//ｙ列（num_row）の列データのforに対応

return 0;
}


本3章の説明へ



9、最小二乗 1/xの-n乗 多項式近似(コード9)
#include <string>
#include <cmath>
#include <fstream>
#include <vector>
#include <algorithm>
#include <functional>
#include <iomanip>
#include <numeric>
using namespace std;

// y=f(x)が低減するデータの近似式、Σ(a/x)のn乗で近似項の係数
int main()
{
// 最小二乗法　多項式近似 データ列数およびデータ行数の制約無し

//0、近似用データの生成
int mm=20;//データ行数
vector<double>xx;
xx.assign(mm,0.0);
vector<double>yy0,yy1,yy2,yy3,yy4;
yy0.assign(mm,0.0);
yy1.assign(mm,0.0);
yy2.assign(mm,0.0);
yy3.assign(mm,0.0);
yy4.assign(mm,0.0);
double xt=0.0;

for(int i=0;i<mm;++i)
{
xt=0.125*i+i*i/500.0-5.7500;
xx[i]=xt;
yy0[i]=2.50+5.0/xt+7.5/xt/xt;
yy1[i]=2.50+5.0/xt+7.5/xt/xt+10.0/xt/xt/xt;
yy2[i]=2.50+5.0/xt+7.5/xt/xt+10.0/xt/xt/xt+12.5/xt/xt/xt/xt+0.15/xt/xt/xt/xt/xt;
yy3[i]=2.50+5.0/xt+7.5/xt/xt+10.0/xt/xt/xt+12.5/xt/xt/xt/xt+0.15/xt/xt/xt/xt/xt+175.0/xt/xt/xt/xt/xt/xt;
yy4[i]=2.50+5.0/xt+7.5/xt/xt+10.0/xt/xt/xt+12.5/xt/xt/xt/xt+0.15/xt/xt/xt/xt/xt+175.0/xt/xt/xt/xt/xt/xt
+500.0/xt/xt/xt/xt/xt/xt/xt/xt;
}
ofstream outxxyy("sample_data.csv",ios::out |ios::binary );//データfile
outxxyy<<"データ名sample兼入力データ作成方法 "<<endl;
outxxyy<<"近似次数最大値"<<endl;
outxxyy<<9<<endl;
outxxyy<<"x値とする列番号の指定(仮設定1…0起点故2列目で1)"<<endl;
outxxyy<<1<<endl;
outxxyy<<"近似対象データの列並び(仮設定2…0起点故3列目以右で2)"<<endl;
outxxyy<<2<<endl;
outxxyy<<"データ番号"<<","<<"x値"<<","<<"y値(以右列)"<<endl;
for(int i=0;i<mm;i++)
outxxyy<<i<<","<<xx[i]<<","<< yy0[i]<<","<<yy1[i]<<","<<yy2[i]<<","<<yy3[i]<<","<<yy4[i]<<endl;
outxxyy.close();

//1、csvファイルから","区切りを利用しデータ読込み //データの読み込み（STL/string使用）
// データ列に制限はない。（上記例のようにxは1列yの横5列であるが、計算可能）
// データ行数に制限はない。
string string0,string1,string2,string3,string4;
string stringdata0;
vector<string>stringdata;
int nmmax;
int x_location;
int y_location;

ifstream input0("sample_data.csv", ios::in | ios::binary);
if(!input0.is_open()) exit(1);
input0>>string0;
input0>>string1;
input0>>nmmax;
input0>>string2;
input0>>x_location;
input0>>string3;
input0>>y_location;
input0>>string4;
while(!input0.eof())
{
input0>>stringdata0;
if(!input0.eof())
{
stringdata.push_back(stringdata0);
}
}
input0.close();
stringdata.resize(stringdata.size());

for(int j=0;j<(int)stringdata.size();j++) stringdata[j]=stringdata[j]+",";
vector<double>xv;
string moji=",";
int kp;
for(int j=0;j<(int)stringdata.size();j++)
{
kp=0;
while(stringdata[j].size())
{
string ::iterator place0=find_first_of(stringdata[j].begin(),stringdata[j].end(),moji.begin(),moji.end());
++kp;
if(place0 !=stringdata[j].end())
{
string s0;
for(int i=0;i<(place0-stringdata[j].begin());i++) s0.insert(s0.end(),stringdata[j][i]);
string s0charcheck=s0;
const char* check0=s0charcheck.c_str();
double xs0;
xs0 = strtod(check0, NULL);
xv.push_back(xs0);
basic_string <char>::iterator strg0_Iter;
strg0_Iter =stringdata[j].erase (stringdata[j].begin(),stringdata[j].begin()+(place0-stringdata[j].begin())+1);
}
else if(place0==(stringdata[j].end()))
{
string s1;
for(int i=0;i<(place0-stringdata[j].begin());i++) s1.insert(s1.end(),stringdata[j][i]);
string s1charcheck=s1;
const char* check1=s1charcheck.c_str();
double xs1;
xs1 = strtod(check1, NULL);
xv.push_back(xs1);
basic_string <char>::iterator strg0_Iter;
strg0_Iter =stringdata[j].erase (stringdata[j].begin(),stringdata[j].begin()+(place0-stringdata[j].begin()));
break;
}
}
}
xv.resize(xv.size());


// 2、xとｙデータ列の組み合わ
// x列は[kp*i+x_location]で変更可
// y列はnum_rowでy_locationからkpまでの数を取る。

vector<double>x,y;
for(int num_row=y_location;num_row<kp;num_row++)//return直前の" ｝"と対応
{

x.clear();
y.clear();

for(int i=0;i<(int)xv.size()/kp;i++)
{
x.push_back(xv[kp*i+x_location]);
y.push_back(xv[kp*i+num_row]);
}
x.resize(x.size());
y.resize(y.size());


//3、 最小二乗法　多項式近似


//3-1、データ数をm、求める多項式の次数をnm
// Σ1/xiの(i<m)の一次、二次、三次、…、2*(nm+1)次までの係数を作成する。最終的には、(nm+1)×(nm+1)正方行列のデータになる。


//2-1、Σ1/xiの2*nm乗のデータ作成
int nm;
int m=(int)xv.size()/kp;
vector<double>v1,v2;

vector<double>v3,v4;
for( nm=2;nm<=nmmax;nm++) //return二つ前の" ｝"と対応
{
v1.clear();
v2.clear();

double mulx;
for(int j=0;j<m;j++)
{
mulx=1.0;
for(int i=0;i<2*nm;i++)
{
mulx*=1.0/x[j];
v1.push_back(mulx);
}
}
v1.resize(v1.size());

double sum_mdata=0.0;
v2.push_back(m);
for(int i=0;i<2*nm;i++)
{
sum_mdata=0.0;
for(int j=0;j<m;j++)
{
sum_mdata+=v1[i+j*2*nm];
}
v2.push_back(sum_mdata);
}
v2.resize(v2.size());

v3.clear();
v4.clear();
double muly=0.0;
for(int j=0;j<m;j++)
{
muly=y[j];
for(int i=0;i<2*nm;i++)
{
muly*=1.0/x[j];
v3.push_back(muly);
}
}
v3.resize(v3.size());

v4.push_back(accumulate(y.begin(),y.end(),0.0));
double resultxy=0.0;
for(int i=0;i<2*nm;i++)
{
resultxy=0.0;
for(int j=0;j<m;j++) resultxy+=v3[i+j*2*nm];
v4.push_back(resultxy);
}
v4.resize(v4.size());

//4、連立方程式の作成と解の算出
vector<double>s,t;
vector<double>ans;
vector<double>ddd;
vector<int>row;
s.clear();
t.clear();
ans.clear();
ddd.clear();
row.clear();

int i,j,k;
int n=nm+1;
t.assign(n,0.0);
for(int j=0;j<n;j++)
{
for(int i=j;i<n+j;i++) s.push_back(v2[i]);
}
s.resize(s.size());
for(int i=0;i<n;i++) t[i]=v4[i];

vector<double>ss,tt;
ss.assign(n*n,0.0);
tt.assign(n,0.0);
ss=s;
tt=t;
ans.assign(n,0.0);
for(i=0;i<n;++i) row.insert(row.end(),1,i);
for(k=0;k<n;k++)
{
ddd.assign(n*n,0.0);
transform(s.begin()+n*k,s.end(),ddd.begin(),fabs);
int inum=max_element(ddd.begin(),ddd.end())-ddd.begin();
int maxvline=(int)((n*k+inum)/n);
int maxvrow=n*k+inum- maxvline*n;
for(i=0;i<n;i++) swap(s[n*k+i],s[n*(maxvline)+i]);
for(i=0;i<n;i++) swap(s[i*n+k],s[i*n+maxvrow]);
swap(t[k],t[maxvline]);
swap(row[k],row[maxvrow]);
for (i=k+1; i<n; ++i)
{
t[i]-=t[k]*s[i*n+k]/s[k*n+k];
for(j=n-1;j>=k;j--)
{
if(fabs(s[k*n+k])>1.0e-15) s[i*n+j] -=s[i*n+k]*s[k*n+j]/s[k*n+k];
}
}
}
double anssum;
anssum=0.0;
for(i=n-1;i>=0;--i)
{
anssum=0.0;
for(j=0;j<n-i-1;j++) anssum=s[n*i+n-1-j]*ans[row[n-1-j]]+anssum;
ans[row[i]]=(t[i]-anssum)/s[n*i+i];
}

//5、計算結果検証用データ作成
vector<double>recal;
double sum;
for(j=0;j<m;j++)
{
sum=0.0;
for(i=1;i<n;i++)
{
sum=sum+ans[i]*pow(1.0/x[j],i);
}
recal.push_back(sum+ans[0]);
}
recal.resize(recal.size());

//6、相関係数R^2の算出
// 1)ε残差とｙ残差による相関係数算出
double average_x=accumulate(x.begin(),x.end(),0.0)/(x.size());
double average_y=accumulate(y.begin(),y.end(),0.0)/(y.size());
double sum_resty=0.0;
for(int i=0;i<(int)y.size();i++) sum_resty+=(y[i]-recal[i])*(y[i]-recal[i]);
double sum_yy=0.0;
for(int i=0;i<(int)y.size();i++) sum_yy+=(y[i]-average_y)*(y[i]-average_y);
double rsqrest_pmcc=1-sum_resty/sum_yy;

ofstream out1("calculation_result.csv",ios::out |ios::binary |ios::app);
out1<<" "<<endl;
out1<<string0<<endl;
out1<<"指定した計算の最高次数(入力値)= "<<nmmax<<endl;
out1<<"多項式　x値の列番号= "<<x_location<<" 近似対象データ列番号= "<<num_row<<endl;
out1<<"多項式適用次数= "<<nm<<endl;
out1<<" "<<endl;
out1<<"計算結果 近似係数の順＝定数、1次…n次　"<<endl;
out1<<"次数"<<","<<"近似項の係数"<<endl;
for(i=0;i<(int)ans.size();i++)
out1<<setw(2)<<i<<","<<setw(15)<<setprecision(10)<<ans[i]<<endl;
out1<<" "<<endl;
out1<<"近似式の検証　"<<endl;
out1<<" i "<<","<<"x値"<<","<<" 再計算値（近似式）"<<","<<" 原値(y)"<<","<<"再計算値-原値"<<endl;
for(i=0;i<(int)recal.size();i++)
out1<<i<<","<<x[i]<<","<<recal[i]<<","<< y[i]<<","<<recal[i]-y[i]<<endl;
out1<<" "<<endl;
out1<<"相関係数R2 "<<","<<rsqrest_pmcc<<endl;
out1<<" "<<endl;
out1<<" "<<endl;
out1.close();

}
}
return 0;

}


本3章の説明へ



10、フーリエ級数展開(コード10)
#include <string>
#include <cmath>
#include <vector>
#include <fstream>
#include <clocale>
#include <numeric>
using namespace std;
#define pi 3.14159265358979

int main(void)
{
setlocale( LC_ALL, "japanese" );

// Fourier級数用入力データ作成
int div_num=223;
ofstream out0("inputdata.csv", ios::out | ios::binary );
double phase1=-pi/4.0;
for(int i=0;i<div_num;i++)
out0<<8.855+50.0*cos(i*2.0*pi/(double)div_num)+100.00*cos(i*2.0*2.0*pi/(double)div_num+phase1)+
20.0*cos(i*4.0*2.0*pi/(double)div_num)+330.0*cos(i*19.0*2.0*pi/(double)div_num-phase1*2.0)<<endl;
out0.close();

// 周期の長さを2dとする。
ifstream input0("inputdata.csv", ios::in | ios::binary);//
if(!input0.is_open()) exit(1);

string stringdata0;
vector<double>v1;//このv1が主役。
double indata;

while(!input0.eof())
{
input0>>stringdata0;
if(!input0.eof())
{
indata=strtod((const char*)stringdata0.c_str(),NULL);
v1.push_back(indata);
}
}
input0.close();
v1.resize(v1.size());


double d=(double)v1.size(); //v1.size()はintのために、double化。
double a0=2.0/d*accumulate(v1.begin(),v1.end(),0.0);//初期値0.0で、v1を加算する。

double c1,c2,fai;
vector<double>cosfourier;
vector<double>sinfourier;
vector<double>phase;

for(int m=0;m<(int)v1.size()/2;m++)//mは次数。pはデータの位置、また次数の範囲は(int)v1.size()/2
{
c1=0.0;
for(int p=0;p<(int)v1.size();p++) c1=c1+2.0/d*v1[p]*sin((double)m*2.0*pi/d*(double)p);
sinfourier.push_back(c1);
c2=0.0;
for(int p=0;p<(int)v1.size();p++) c2=c2+2.0/d*v1[p]*cos((double)m*2.0*pi/d*(double)p);
cosfourier.push_back(c2);
fai=180.0/pi*atan(c1/(c2+1.0e-10));
phase.push_back(fai);
}
cosfourier.resize(cosfourier.size());
sinfourier.resize(sinfourier.size());
phase.resize(phase.size());

vector<double>cosy,siny,sum1;
for(int p=0;p<(signed)v1.size();p++)
{
c1=0.0;
c2=0.0;
for(int m=1;m<(int)cosfourier.size();m++) c2=c2+cosfourier[m]*cos((double)m*2.0*pi/d*(double)p);
cosy.push_back(c2);

for(int m=1;m<(int)sinfourier.size();m++) c1=c1+sinfourier[m]*sin((double)m*2.0*pi/d*(double)p);
siny.push_back(c1);
sum1.push_back(c1+c2+a0/2);
}
cosy.resize(cosy.size());
siny.resize(siny.size());
sum1.resize(sum1.size());


ofstream out2c("outfourier.csv", ios::out | ios::binary );
out2c<<" "<<endl;
out2c<<" a0/2= "<<a0/2<<endl;
for(int i=0;i<(int)cosfourier.size();i++) out2c<<" cos 次数["<<i<<"]= "<<","<<cosfourier[i]<<endl;
out2c<<" "<<endl;
for(int i=0;i<(int)sinfourier.size();i++) out2c<<" sin 次数["<<i<<"]= "<<","<<sinfourier[i]<<endl;
out2c<<" "<<endl;
for(int i=0;i<(int)phase.size();i++) out2c<<" 位相差 次数["<<i<<"]= "<<","<<phase[i]<<endl;
out2c<<" "<<endl;
out2c<<" p "<<","<<"�把os[i]"<<","<<" v1[i] "<<endl;
for(int i=0;i<(int)cosy.size();i++) out2c<<i<<","<<cosy[i]<<","<<v1[i]<<endl;
out2c<<" "<<endl;
out2c<<" p "<<","<<"�敗in[i]"<<","<<" v1[i] "<<endl;
for(int i=0;i<(int)siny.size();i++) out2c<<i<<","<<siny[i]<<","<<v1[i]<<endl;
out2c<<" "<<endl;
out2c<<" p "<<","<<"a0/2+�把os+�敗in"<<","<<" v1[i] "<<endl;
for(int i=0;i<(int)siny.size();i++) out2c<<i<<","<<sum1[i]<<","<<v1[i]<<endl;
out2c.close();
return 0;
}


本3章の説明へ

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　copyright©2009-20xx Nakadaka&Kobaru all rights reserved.