---

_________________________＿＿＿＿__＿_____＿＿＿________計算コードの具体例______________________________________________
1、珪素鋼板 BH 特性に対する補間計算(コード 1 )
#include <cmath>
#include <iostream>
#include <fstream>
#include <vector>
#include <algorithm>
#include <functional>
#include <iomanip>
#include<numeric>
using namespace std;

const double pi=3.14159265358979;

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

double b30(int);
double h30(int);


int main()
{

int jbh;
vector<double>bhdata;
vector<double>bdata;
vector<double>hdata;

for(jbh=0;jbh<79;jbh++)
{
bdata.push_back(b30(jbh));
hdata.push_back(h30(jbh));
}

hdata.resize(hdata.size());
bdata.resize(bdata.size());

//最初に、hを入れて、次にhの先頭個所に、bのデータを入れる。
vector<double>::iterator p1=bhdata.begin();
vector<double>::iterator outstart=hdata.begin();
vector<double>::iterator outend=hdata.end();
bhdata.insert(p1,outstart,outend);

//hの先頭個所に、bのデータを入れる。
vector<double>::iterator p2=bhdata.begin();
vector<double>::iterator out2start=bdata.begin();
vector<double>::iterator out2end=bdata.end();
bhdata.insert(p2,out2start,out2end);

bhdata.resize(bhdata.size());
// 結局　bhdataには、bが始めの部分にあり、次の部分にhがある。



// 補間行うx値
vector<double>v;
int num=10;
v.assign(num,0.0);
v[0]=50.0;
v[1]=800.0;
v[2]=8000.0;
v[3]=12000.0;
v[4]=13000.0;
v[5]=14500.0;
v[6]=16500.0;
v[7]=18000.0;
v[8]=22000.0;
v[9]=26000.0;

// 本計算のための前準備
int m=(int)bdata.size();
vector<double>x(m,0.0);
x=bdata;
vector<double>y(m,0.0);
y=hdata;
double x0;//補間点


// Lagrange補間
for(int i=0;i<num;i++)//補間点データの繰り返し
{
x0=v[i];
vector<double>lag1(m,0.0);
vector<double>lag2(m,0.0);
double product0,product1;
vector<double>productx(m,0.0);

for(int j=0;j<m;j++)
{
for(int i=0;i<m;i++)
{
lag1[i]=x[j]-x[i];
lag2[i]=x0-x[i];
}
replace(lag1.begin(),lag1.end(),0.0,1.0);
product0=accumulate(lag1.begin(),lag1.end(),1.0,multiplies<double>());
replace(lag2.begin(),lag2.end(),x0-x[j],1.0);
product1=accumulate(lag2.begin(),lag2.end(),1.0,multiplies<double>());
productx[j]=product1/product0;
}
double lagrange_value=inner_product(y.begin(),y.end(),productx.begin(),0.0);

ofstream out2("out_lagrange.csv",ios::out |ios::binary |ios::app);//app追加
out2<<" Lagarange interaption "<<endl;
out2<<" y_interaption value= "<<lagrange_value<<" at x= "<<x0<<endl;
out2<<" "<<endl;
out2.close();
}


// Lagrange補間の方程式の求根形
vector<double>yyy(m,0.0);// add for this caluculation
yyy=hdata;

double eps=1.0e-6;
vector<double>lag1(m,0.0);
vector<double>lag2(m,0.0);
double product0,product1;
vector<double>producty(m,0.0);
vector<double>xxx(m,0.0);

for(int i=0;i<num;i++)// 補間点データの繰り返し
{
x0=v[i];
for(int i=0;i<m;++i) xxx[i]=x[i]-x0;
vector<double>v0(m,0.0);
transform(xxx.begin(),xxx.end()-1,xxx.begin()+1,v0.begin(),multiplies<double>());
int near0=find_if(v0.begin(),v0.end(),minusvalue)-v0.begin();//x0に近いx[i]の位置算出。
double a=yyy[near0];//yを計算する区間算出
double b=yyy[near0+1];//yを計算する区間算出
int n=(int)(log(fabs(b-a)/eps)/log(2.0)+0.5);//繰り返し回数で区間を分割する事にする。
double h=(b-a)/n;
vector<double>lagrange_y(n+1,0.0);
double ypitch;

int near01=0;
int jk=0;
int ik=0;
for(int kk=0;kk<200;++kk)//収束数値の精度upのために200とした。
{
h=(b-a)/n;
jk++;
for(int k=0;k<=n;k++)
{
ypitch=a+k*h;
ik++;
for(int j=0;j<m;j++)
{
for(int i=0;i<m;i++)
{
lag1[i]=yyy[j]-yyy[i];
lag2[i]=ypitch-yyy[i];
}
replace(lag1.begin(),lag1.end(),0.0,1.0);
product0=accumulate(lag1.begin(),lag1.end(),1.0,multiplies<double>());
replace(lag2.begin(),lag2.end(),ypitch-yyy[j],1.0);
product1=accumulate(lag2.begin(),lag2.end(),1.0,multiplies<double>());
producty[j]=product1/product0;
}
lagrange_y[k]=inner_product(x.begin(),x.end(),producty.begin(),0.0)-x0;

if(fabs(lagrange_y[k])<eps*0.01)//x0がベースのx[i]と同じ値の時に適用となる。最初にゼロに収束。
{
ofstream out3("out_lagarange_equation1.csv",ios::out |ios::binary |ios::app);
out3<<" 細分化step jk= "<<jk<<" step ik= "<<ik<<endl;
out3.precision(12);
out3<<" y_interaption value= "<<ypitch<<" at x= "<<x0<<endl;
out3<<" Note: y_interaption_value may be in inputed y_data. please, check x&y_data "<<endl;
out3.close();
break;
}
}
vector<double>fabslag(n+1,0.0);
for(int i=0;i<=n;++i) fabslag[i]=fabs(lagrange_y[i]);
int near01=min_element(fabslag.begin(),fabslag.end())- fabslag.begin();
fabslag.clear();

b=a+(near01+1)*h;
a=a+(near01)*h;
if(fabs(lagrange_y[near01])<eps) break;
}
ofstream out2("out_lagarange_equation.csv",ios::out |ios::binary |ios::app);
out2<<" Lagarange_equation interaption "<<endl;
out2<<" 細分化step jk= "<<jk<<" step ik= "<<ik<<endl;
out2.precision(12);
out2<<" y_interaption value= "<<ypitch<<" at x= "<<x0<<endl;
out2<<" "<<endl;
out2.close();
}


// Aitken補間
vector<double> xx(m,0.0);
vector<double> yy(m,0.0);
xx=bdata;
yy=hdata;
vector<double>xt(m,0.0),yt(m,0.0);
xt=xx;
yt=yy;

for(int i=0;i<num;i++)
{
x0=v[i];
vector<double>xat(m,-x0);
vector<double>xxt(m,0.0);
transform(xt.begin(),xt.end(),xat.begin(),xxt.begin(),plus<double>());//補間点xとx[i]の差
vector<double>xxxt(m,0.0);
transform(xxt.begin(),xxt.end(),xxxt.begin(),fabs);//±の差を絶対値化(xxxt)し、最小順に並べる準備

vector<double>t1(m,0.0);
t1=xxxt;// 絶対値のxxxtを最小順に並べ変えるので、t1に保存
sort(xxxt.begin(),xxxt.end(),less<double>());//less<double>()により最小順に並変え。

//隣合わせが同じ数値の場合、同じ標本点のデータが続けてx,yに挿入され、計算エラーが生じる。
//このために、僅か(1.0000001倍)に補間点を動かして、補間値を求める。
vector<double>t2(m,0.0);
adjacent_difference(xxxt.begin(),xxxt.end(),t2.begin());// 同じ数値が隣合わせにないか調査
int number0=count_if(t2.begin(),t2.end(),bind2nd(equal_to<double>(),0.0));
if(number0 !=NULL)
{
xxxt.assign(m,0.0);
t1.assign(m,0.0);
x0=x0*(1.0+0.000001);
vector<double>xat(m,-x0);
vector<double>xxt(m,0.0);
transform(xt.begin(),xt.end(),xat.begin(),xxt.begin(),plus<double>());
transform(xxt.begin(),xxt.end(),xxxt.begin(),fabs);
t1=xxxt;
sort(xxxt.begin(),xxxt.end(),less<double>());
vector<double>t2(m,0.0);
adjacent_difference(xxxt.begin(),xxxt.end(),t2.begin());
int number01=count_if(t2.begin(),t2.end(),bind2nd(equal_to<double>(),0.0));
if(number01 !=NULL)
{
ofstream out01(" alarm1.csv",ios::out |ios::binary |ios::app );
out01<<" Please, sample data are always same pitch "<<endl;
out01<<" original data x0 = "<<x0/(1.0-0.000001)<<endl;
out01<<" same pitch count number(original condition) "<<number0<<endl;
out01<<" 1.000001 times changed data x0 = "<<x0<<endl;
out01<<" same pitch count number(changed condition) "<<number01<<endl;
out01<<" data before_sort 　ti "<<endl;
for(int i=0;i<m;i++) out01<<i<<","<<t1[i]<<endl;
out01.close();
break;
}
}
vector<double>x,y;//補間計算用の新データvector
for(int i=0;i<m;++i)
{
int position1=find(t1.begin(),t1.end(),xxxt[i])- t1.begin();// 並び変え前の最小値の位置の算出
x.push_back(xx[position1]);
y.push_back(yy[position1]);
}
x.resize(x.size());
y.resize(y.size());

//4、 aiken補間の計算　次数をi、標本点をjで操作
vector<double>pkn(m*m,0.0);
double eps=1.0e-3;
int kk;
for(int j=0;j<m;j++) pkn[j]=y[j];
for(int i=1;i<m;i++)
{
for(int j=i;j<m;j++)
{
pkn[i*m+j]=(pkn[(i-1)*m+j]*(x0-x[(i-1)])-pkn[(i-1)*m+(i-1)]*(x0-x[j]))/(x[j]-x[(i-1)]);

}
if(fabs((pkn[(i-1)*m+(i-1)]-pkn[i*m+i])/pkn[i*m+i])<eps) break;
kk=i;
}
ofstream out0("out_aitken.csv",ios::out |ios::binary|ios::app );
out0<<" Aitken interaption "<<endl;
out0<<" interpolation order = "<<kk<<endl;
out0<<" y_interaption value= "<<pkn[kk*m+kk]<<" at x= "<<x0<<endl;
out0<<" "<<endl;
out0.close();
}

//spline 補間
vector<double>u;
u.assign(bhdata.size(),0.0);
u=bhdata;

int j,i;
int spn;
double dx;
spn=u.size()/2;
vector<double>w;
w.assign(v.size(),0.0);
vector<double>s0;
s0.assign(4*spn,0.0);
for(i=0;i<2*spn;i++) s0[i]=u[i];
adjacent_difference(s0.begin(),s0.begin()+spn,s0.begin()+2*spn);
adjacent_difference(s0.begin()+spn,s0.begin()+2*spn,s0.begin()+3*spn);
vector<double>tsp,ssp;
tsp.assign(spn-2,0.0);
for(i=0;i<spn-2;i++) tsp[i]=s0[3*spn+2+i]/s0[2*spn+2+i]-s0[3*spn+1+i]/s0[2*spn+1+i];
ssp.assign((spn-2)*(spn-2),0.0);
for(i=0;i<spn-2;i++) ssp[(spn-2)*i+i]=(s0[2*spn+2+i]+s0[2*spn+1+i])/3.0;
for(i=0;i<spn-3;i++) ssp[(spn-2)*i+i+1]=s0[2*spn+2+i]/6.0;
for(i=0;i<spn-3;i++) ssp[(spn-2)*(i+1)+i]=s0[2*spn+2+i]/6.0;

vector<double>t,s;
s=ssp;
t=tsp;
int n=spn-2;
vector<double>ss,tt;
ss.assign(n*n,0.0);
tt.assign(n,0.0);
ss=s;
tt=t;
int k;
//vector<double>x;
vector<double>ans;
ans.assign(n,0.0);
vector<double>ddd;
vector<int>row;
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=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];
}
for(i=0;i<(int)ans.size();i++) cout<<setw(2)<<i<<","<<setw(15)<<setprecision(10)<<ans[i]<<endl;

ans.insert(ans.begin(),1,0.0);
ans.insert(ans.end(),1,0.0);

int mnum,ii;
vector<double>::iterator dxgreater;
for(ii=0;ii<(int)v.size();ii++)
{
dx=v[ii];
dxgreater=search_n(s0.begin(),s0.begin()+spn,1,dx,greater<double>());
for(i=0;i<=(int)s0.size()/4;i++) if(s0[i]- *dxgreater ==0.00) mnum=i;
if(mnum>=(int)s0.size()/4) mnum=s0.size()/4-1;
double u0=s0[mnum]-dx;
double u1=dx-s0[mnum-1];
double v0=ans[mnum-1]*(pow(u0,3)/s0[2*spn+mnum]-u0*s0[2*spn+mnum])/6.0;
double v1=ans[mnum]*(pow(u1,3)/s0[2*spn+mnum]-u1*s0[2*spn+mnum])/6.0;
double w0=v0+v1+(u0*s0[spn+mnum-1]+u1*s0[spn+mnum])/s0[2*spn+mnum];
w[ii]=w0;
}

ofstream cout0("out_spline.csv",ios::out |ios::binary);
cout0<<" spline interaption "<<endl;
cout0<<"補間点番号 i"<<","<<" 補間値w[i]"<<","<<" at 補間点v[i]"<<endl;
for(i=0;i<(int)w.size();i++) cout0<<i<<","<<w[i]<<","<<v[i]<<endl;
cout0.close();

return 0;

}

//珪素鋼板のB-H特性
double b30(int j)
{
double b30mm[79]={ 0.0,10.80,19.45,21.61,64.84,108.08,151.31,216,376,585,864,1247,1767,2450,3278,4188,6023,7598,8842,9794,
10524,11518,12184,12648,12990,13252,13465,13636,13784,13907,14018,14237,14404,14663,14854,15007,15135,
15246,15345,16009,16449,16795,17094,17808,18180,
18300,18500,18820,19000,19400,19500,19840,20000,20250,20600,20920,21200,21480,21980,22350,22800,23200,
23400,23650,24000,24500,25200,25900,26500,27250,29400,34200,38100,47990,47995,47997,47998,47999,48000};

return b30mm[j];
}

double h30(int j)
{
double h30mm[79]={ 0.0,0.005,0.009,0.01,0.03,0.05,0.07,0.099925,0.14988,0.19981,0.24975,0.2997,0.34965,0.39929,0.44914,
0.49889,0.5985,0.69823,0.79795,0.89774,0.99752,1.193,1.3926,1.5923,1.7923,1.9919,2.1922,2.3919,2.5917,
2.7916,2.9915,3.4914,3.9782,4.9777,5.9768,6.977,7.9762,8.9752,9.974,19.92,29.918,39.775,49.764,79.749,
99.732,110.001,125.001,150.001,163.001,200.001,210.001,250.001,265.01,300.01,350.001,400.01,450.001,
500.001,600.01,700,800,900,950,1000,1100,1200.12,1350.13,1500.15,1650.16,1800.1,2292.3,3369.2,
4138.4,6443.668066,6444.908043,6445.405426,6445.653617,6445.901809,6446.15};

return h30mm[j];
}

本章の説明へ


2、直交関数 Legendre 多項式による補間計算 (コード 2 )
#include <cstdio>
#include <cmath>
#include <fstream>
#include <vector>
#include <iomanip>
#include <numeric>
#include <iostream>
using namespace std;

vector<double> Legendre_function();
double b30(int);
double h30(int);


int main()
{

double xj=26000.0 ;//補間

int div0;
int div;
int ak;
int parameter;
int first; //区分の最初の番号

if(b30(0)<=xj && xj<b30(4)) first=0, div0=6, ak=2, parameter=3,div=div0/parameter;//0-4
if(b30(4)<=xj && xj<b30(8)) first=1, div0=6, ak=2, parameter=3,div=div0/parameter;//4-8
if(b30(8)<=xj && xj<b30(9)) first=7, div0=4, ak=2, parameter=4,div=div0/parameter;//8-9
if(b30(9)<=xj && xj<b30(10)) first=1, div0=18,ak=7, parameter=6,div=div0/parameter;//9-10
if(b30(10)<=xj && xj<b30(11)) first=7, div0=4, ak=2, parameter=4,div=div0/parameter;//10-11
if(b30(11)<=xj && xj<b30(12)) first=9, div0=10,ak=5, parameter=2,div=div0/parameter;//11-12
if(b30(12)<=xj && xj<b30(19)) first=12,div0=8, ak=5, parameter=8,div=div0/parameter;//12-19
if(b30(19)<=xj && xj<b30(20)) first=13,div0=8, ak=5, parameter=8,div=div0/parameter;//19-20
if(b30(20)<=xj && xj<b30(22)) first=19,div0=4, ak=3, parameter=4,div=div0/parameter;//20-22
if(b30(22)<=xj && xj<b30(23)) first=19,div0=12,ak=5, parameter=4,div=div0/parameter;//22-23
if(b30(23)<=xj && xj<b30(25)) first=23,div0=6, ak=3, parameter=6,div=div0/parameter;//23-25
if(b30(25)<=xj && xj<b30(37)) first=25,div0=12,ak=4, parameter=4,div=div0/parameter;//25-37
if(b30(37)<=xj && xj<b30(42)) first=34,div0=8, ak=4, parameter=8,div=div0/parameter;//37-42
if(b30(42)<=xj && xj<b30(44)) first=38,div0=6, ak=3, parameter=6,div=div0/parameter;//42-44
if(b30(44)<=xj && xj<b30(45)) first=44,div0=6, ak=3, parameter=6,div=div0/parameter;//44-45
if(b30(45)<=xj && xj<b30(46)) first=45,div0=8, ak=6, parameter=8,div=div0/parameter;//45-46
if(b30(46)<=xj && xj<b30(48)) first=45,div0=3, ak=2, parameter=3,div=div0/parameter;//46-48
if(b30(48)<=xj && xj<b30(49)) first=47,div0=4, ak=2, parameter=4,div=div0/parameter;//48-49
if(b30(49)<=xj && xj<b30(50)) first=49,div0=8, ak=2, parameter=8,div=div0/parameter;//48-50
if(b30(50)<=xj && xj<b30(51)) first=47,div0=4, ak=2, parameter=4,div=div0/parameter;//50-51
if(b30(51)<=xj && xj<b30(52)) first=49,div0=3, ak=2, parameter=3,div=div0/parameter;//51-52
if(b30(52)<=xj && xj<b30(53)) first=37,div0=16,ak=9, parameter=8,div=div0/parameter;//52-53
if(b30(53)<=xj && xj<b30(55)) first=52,div0=5, ak=3, parameter=5,div=div0/parameter;//53-55
if(b30(55)<=xj && xj<b30(57)) first=55,div0=2, ak=2, parameter=2,div=div0/parameter;//55-57
if(b30(57)<=xj && xj<b30(58)) first=46,div0=18,ak=9, parameter=9,div=div0/parameter;//57-58
if(b30(58)<=xj && xj<b30(60)) first=56,div0=3, ak=2, parameter=3,div=div0/parameter;//58-60
if(b30(60)<=xj && xj<b30(63)) first=55,div0=8, ak=6, parameter=8,div=div0/parameter;//60-63
if(b30(63)<=xj && xj<b30(65)) first=56,div0=8, ak=4, parameter=8,div=div0/parameter;//63-65
if(b30(65)<=xj && xj<b30(68)) first=62,div0=6, ak=3, parameter=6,div=div0/parameter;//65-68
if(b30(68)<=xj && xj<b30(69)) first=59,div0=10,ak=4, parameter=5,div=div0/parameter;//68-69
if(b30(69)<=xj && xj<b30(70)) first=68,div0=4, ak=3, parameter=4,div=div0/parameter;//69-70
if(b30(70)<=xj && xj<b30(73)) first=69,div0=3, ak=2, parameter=3,div=div0/parameter;//70-73
if(b30(73)<=xj && xj<b30(78)) first=73,div0=5, ak=4, parameter=5,div=div0/parameter;//73-78


int last=first+div0; //区分の最後の番号

double a=b30(first);
double b=b30(last);

double x=0.0;
double bj=0.0;

vector<double>g_x;
for(int j=0;j<first;j++) g_x.push_back(0.0);

for(int j=first;j<last+1;j++)
{
bj=(2.0*b30(j)-(b+a))/(b-a);
g_x.push_back(bj);
}
g_x.resize(g_x.size());


vector<double>g_value;
for(int j=0;j<80;j++) g_value.push_back(h30(j));
g_value.resize(g_value.size());


vector<double>g_f;
g_f=g_value;
g_f.resize(g_f.size());


vector<double>Interpolation_an(ak+1,0.0);

vector<double>hx;
hx=Legendre_function();
hx.resize(hx.size());

vector<double>fx;

for(int L=ak;L>=0;L--)
{
fx.clear();
int m=L;
for(int i=L*(L+1)/2;i<m+L*(L+1)/2+1;i++) fx.push_back(hx[i]);
fx.resize(fx.size());

double factor=(2.0*L+1)/2.0;

double f_xn00=0.0;//分割始点
double f_xn10=0.0;//幅の1/10点
double f_xn11=0.0;//幅の1/9点
double f_xn13=0.0;//幅の1/8点
double f_xn14=0.0;//幅の1/7点
double f_xn17=0.0;//幅の1/6点
double f_xn20=0.0;//幅の1/5 2/10点
double f_xn22=0.0;//幅の2/9点
double f_xn25=0.0;//幅の1/4 2/8点
double f_xn29=0.0;//幅の2/7点
double f_xn30=0.0;//幅の3/10点
double f_xn33=0.0;//幅の1/3 2/6 3/9点
double f_xn38=0.0;//幅の3/8点
double f_xn40=0.0;//幅の2/5 4/10点
double f_xn43=0.0;//幅の3/7点
double f_xn44=0.0;//幅の4/9点
double f_xn50=0.0;//中間点 1/2 2/4 3/6 4/8 5/10
double f_xn56=0.0;//幅の5/9点
double f_xn57=0.0;//幅の4/7点
double f_xn60=0.0;//幅の3/5 6/10点
double f_xn63=0.0;//幅の5/8点
double f_xn66=0.0;//幅の2/3 4/6 6/9点
double f_xn70=0.0;//幅の7/10点
double f_xn71=0.0;//幅の5/7点
double f_xn75=0.0;//幅の3/4 6/8点
double f_xn78=0.0;//幅の7/9点
double f_xn80=0.0;//幅の4/5 8/10点
double f_xn83=0.0;//幅の5/6点
double f_xn86=0.0;//幅の6/7点
double f_xn88=0.0;//幅の7/8点
double f_xn89=0.0;//幅の8/9点
double f_xn90=0.0;//幅の9/10点
double f_xned=0.0;//終点

double g_xn00=0.0;//分割始点
double g_xn10=0.0;//幅の1/10点
double g_xn11=0.0;//幅の1/9点
double g_xn13=0.0;//幅の1/8点
double g_xn14=0.0;//幅の1/7点
double g_xn17=0.0;//幅の1/6点
double g_xn20=0.0;//幅の1/5 2/10点
double g_xn22=0.0;//幅の2/9点
double g_xn25=0.0;//幅の1/4 2/8点
double g_xn29=0.0;//幅の2/7点
double g_xn30=0.0;//幅の3/10点
double g_xn33=0.0;//幅の1/3 2/6 3/9点
double g_xn38=0.0;//幅の3/8点
double g_xn40=0.0;//幅の2/5 4/10点
double g_xn43=0.0;//幅の3/7点
double g_xn44=0.0;//幅の4/9点
double g_xn50=0.0;//中間点 1/2 2/4 3/6 4/8 5/10
double g_xn56=0.0;//幅の5/9点
double g_xn57=0.0;//幅の4/7点
double g_xn60=0.0;//幅の3/5 6/10点
double g_xn63=0.0;//幅の5/8点
double g_xn66=0.0;//幅の2/3 4/6 6/9点
double g_xn70=0.0;//幅の7/10点
double g_xn71=0.0;//幅の5/7点
double g_xn75=0.0;//幅の3/4 6/8点
double g_xn78=0.0;//幅の7/9点
double g_xn80=0.0;//幅の4/5 8/10点
double g_xn83=0.0;//幅の5/6点
double g_xn86=0.0;//幅の6/7点
double g_xn88=0.0;//幅の7/8点
double g_xn89=0.0;//幅の8/9点
double g_xn90=0.0;//幅の9/10点
double g_xned=0.0;//終点


double Interpolation_ort=0.0;


for(int j=0;j<div;j++)
{

vector<double>p_xn00(m+1,0.0);//分割始点
vector<double>p_xn10(m+1,0.0);//幅の1/10点
vector<double>p_xn11(m+1,0.0);//幅の1/9点
vector<double>p_xn13(m+1,0.0);//幅の1/8点
vector<double>p_xn14(m+1,0.0);//幅の1/7点
vector<double>p_xn17(m+1,0.0);//幅の1/6点
vector<double>p_xn20(m+1,0.0);//幅の1/5 2/10点
vector<double>p_xn22(m+1,0.0);//幅の2/9点
vector<double>p_xn25(m+1,0.0);//幅の1/4 2/8点
vector<double>p_xn29(m+1,0.0);//幅の2/7点
vector<double>p_xn30(m+1,0.0);//幅の3/10点
vector<double>p_xn33(m+1,0.0);//幅の1/3 2/6 3/9点
vector<double>p_xn38(m+1,0.0);//幅の3/8点
vector<double>p_xn40(m+1,0.0);//幅の2/5 4/10点
vector<double>p_xn43(m+1,0.0);//幅の3/7点
vector<double>p_xn44(m+1,0.0);//幅の4/9点
vector<double>p_xn50(m+1,0.0);//中間点 1/2 2/4 3/6 4/8 5/10
vector<double>p_xn56(m+1,0.0);//幅の5/9点
vector<double>p_xn57(m+1,0.0);//幅の4/7点
vector<double>p_xn60(m+1,0.0);//幅の3/5 6/10点
vector<double>p_xn63(m+1,0.0);//幅の5/8点
vector<double>p_xn66(m+1,0.0);//幅の2/3 4/6 6/9点
vector<double>p_xn70(m+1,0.0);//幅の7/10点
vector<double>p_xn71(m+1,0.0);//幅の5/7点
vector<double>p_xn75(m+1,0.0);//幅の3/4 6/8点
vector<double>p_xn78(m+1,0.0);//幅の7/9点
vector<double>p_xn80(m+1,0.0);//幅の4/5 8/10点
vector<double>p_xn83(m+1,0.0);//幅の5/6点
vector<double>p_xn86(m+1,0.0);//幅の6/7点
vector<double>p_xn88(m+1,0.0);//幅の7/8点
vector<double>p_xn89(m+1,0.0);//幅の8/9点
vector<double>p_xn90(m+1,0.0);//幅の9/10点
vector<double>p_xned(m+1,0.0);//終点

if(parameter==1)
{
// 台形則
for(int i=0;i<m+1;i++)
{
p_xn00[i]=pow(g_x[j+0+first],i);
p_xned[i]=pow(g_x[j+1+first],i);
}

f_xn00=inner_product(fx.begin(),fx.end(),p_xn00.begin(),0.0, plus<double>(),multiplies<double>());
f_xned=inner_product(fx.begin(),fx.end(),p_xned.begin(),0.0, plus<double>(),multiplies<double>());

g_xn00=g_f[1*j+0+first];
g_xned=g_f[1*j+1+first];

Interpolation_ort=Interpolation_ort+(f_xn00*g_xn00+f_xned*g_xned)*(g_x[j+1+first]-g_x[j+0+first])/2.0;
}


if(parameter==2)
{
//simpson1/3則
for(int i=0;i<m+1;i++)
{
p_xn00[i]=pow(g_x[2*j+0+first],i);
p_xn50[i]=pow(g_x[2*j+1+first],i);
p_xned[i]=pow(g_x[2*j+2+first],i);
}

f_xn00=inner_product(fx.begin(),fx.end(),p_xn00.begin(),0.0, plus<double>(),multiplies<double>());
f_xn50=inner_product(fx.begin(),fx.end(),p_xn50.begin(),0.0, plus<double>(),multiplies<double>());
f_xned=inner_product(fx.begin(),fx.end(),p_xned.begin(),0.0, plus<double>(),multiplies<double>());

g_xn00=g_f[2*j+0+first];
g_xn50=g_f[2*j+1+first];
g_xned=g_f[2*j+2+first];

double x0=g_x[2*j+0+first];
double x1=g_x[2*j+1+first];
double x2=g_x[2*j+2+first];

Interpolation_ort=Interpolation_ort+(x2-x0)/(x1-x0)*(3.0*(x1-x0)-(x2-x0))/6.0*f_xn00*g_xn00
+(x2-x0)*(x2-x0)*(x2-x0)/(x2-x1)/(x1-x0)/6.0*f_xn50*g_xn50
+(x2-x0)/(x2-x1)*(3.0*(x2-x1)-(x2-x0))/6.0*f_xned*g_xned;
}


// 以下　第9章 数値計算7 数値積分(離散系)に　記載するので　あるいは、後日記載するので、削除




// 以上　第9章 数値計算7 数値積分(離散系)に　記載するので　あるいは、後日記載するので、削除


}



}//div

Interpolation_an[L]=Interpolation_ort*factor;

}//fx

fx.clear();


int jk=div0+1;




// interpolation
double Interpolation_value=0.0;


for(int L=ak;L>=0;L--)
{
fx.clear();
int m=L;
for(int i=L*(L+1)/2;i<m+L*(L+1)/2+1;i++) fx.push_back(hx[i]);
fx.resize(fx.size());

x=(2.0*xj-(b+a))/(b-a);

for(int i=0;i<m+1;i++)
{
Interpolation_value=Interpolation_value+Interpolation_an[L]*fx[i]*pow(x,i);

}
}


fstream outinterp("bh_Interpolation.csv",ios::out|ios::binary|ios::app);
outinterp<<"積分域 No."<<","<<"Legendre_Pn"<<","<<"積分 次数"<<","<<"x in [-1_1]"<<","<<"xj"<<","<<"Intrpolation_value"<<endl;
outinterp<<first<<","<<ak<<","<<parameter<<","<<setw(3)<<x<<","<<xj<<","<<setprecision(20)<<Interpolation_value<<endl;
outinterp<<" "<<endl;
outinterp.close();


//}


return 0;
}

vector<double> Legendre_function()
{
//Legendre関数p　
vector<long double >p0(0+1,0.0);
p0[0]=1.0;

vector<long double >p1(1+1,0.0);
p1[0]=0.0;
p1[1]=1.0;

// 以下　第7章 数値計算5 直交関数に　記載しているので　あるいは、後日記載するので、削除



// 以上　第7章 数値計算5 直交関数に　記載しているので　あるいは、後日記載するので、削除

p31[31]=14544636039226909/67108864.0;



int k;
vector<double>hx;
k=0;for(int i=0;i<k+1;i++) hx.push_back(p0[i]);
k=1;for(int i=0;i<k+1;i++) hx.push_back(p1[i]);

//　以下　繰り返しのために　削除


//　以上　繰り返しのために　削除
k=31;for(int i=0;i<k+1;i++) hx.push_back(p31[i]);
hx.resize(hx.size());

return hx;

}

double b30(int j)
//　以下　コード1と同じ故に　削除


//　以上　コード1と同じ故に　削除
}


double h30(int j)
//　以下　コード1と同じ故に　削除


//　以上　コード1と同じ故に　削除
}


本章の説明へ

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　copyright©2009-2012 Nakadaka&Kobaru all rights reserved.