//
#include "stdafx.h"
#include "iostream"
#include "fstream"
#include "conio.h"
#include "cmath"
using namespace std;
const int N=2;//Qn(x)中取到的n的值
const int M=N*(N+3)/2;//方程组的阶数
const double PI=3.1415926;
class EPC
{
private:
int i,j,k,l;//循环下标
double Alpha,S0,K22;//系统本身的参数
double eps,t,h,Imax;//eps：精度要求，t：控制变量，h：增量初值，Imax：最大迭代次数
double a[N+1][N+1];//假设的PDF方程的系数未知量的存放矩阵
double x[M];//把a[i][j]矩阵转化而来，用于牛顿下山法中
double f[M];//通过FPK方程形成的，形如Eq.16的方程组。
double EX1[3*N+1],EX2[3*N+1];//X1与X2的K阶矩的存放数组
double Sigma1,Sigma2;
double A[M][M],B[M];
public:
void input1();//通过文档输入系统本身的各种参数
void input2();//通过文档输入假设的PDF方程的系数未知量的存放矩阵的初始值
void transfer1();//把a[i][j]矩阵转化为牛顿下山法中的x矩阵。
void transfer2();//把牛顿下山法中的x矩阵转化为牛顿下山法中的a[i][j]矩阵
double factorial(int k);//计算K的阶乘
void Kthmoment();//计算K阶矩
void func();//通过FPK方程形成的，形如Eq.16的方程组的值。
void netn_root (); //执行拟牛顿法
void gauss (); //全选主元Gauss消去法
void output();
};
void EPC::input1()//通过文档输入系统本身的各种参数
{
cout << "请准备好“Parameters1.txt”,依次存放Alpha,S0,K22,eps,t,h,Imax" << endl;
ifstream fin1("Parameters1.txt");
ofstream fout1("Parameters1-check.txt");
if(!fin1)
{
cout << "不能打开此文件!" << endl;
getch();
exit(1);
}
fin1 >> Alpha >> S0 >> K22 >> eps >> t >> h >> Imax;
//fout1 << Alpha << " " << S0 << " " << K22 << " " << eps << " " << t << " " << h << " " << Imax <<endl;
}
void EPC::input2()//通过文档输入假设的PDF方程的系数未知量的存放矩阵的初始值
{
cout << "请准备好“Parameters2.txt”,依次存放假设的PDF方程的系数未知量的存放矩阵的初始值" << endl;
ifstream fin2("Parameters2.txt");
ofstream fout("Parameters2-check.txt");
if(!fin2)
{
cout << "不能打开此文件!" << endl;
getch();
exit(1);
}
for(i=0;i<=N;i++)
for(j=0;j<=N;j++)
a[i][j]=0;
for(i=1;i<=N;i++)
{
for(j=0;j<=i;j++)
{
fin2 >> a[i][j];
}
}
/*
for(i=0;i<=N;i++)
{
for(j=0;j<=N;j++)
fout << "a[" << i << "][" << j << "]=" << a[i][j] << " ";
fout << endl;
}
*/
}
void EPC::transfer1()//把a[i][j]矩阵转化为牛顿下山法中的x矩阵。
{
double X[(N+2)*(N+2)];
ofstream fout("transfer1-check.txt");
int temp[M+1];
int s=0;
for(i=1;i<=N;i++)
{
for(j=0;j<=i;j++)
{
X[i*(N+2)+j]=a[i][j];
temp[s]=(i*(N+2)+j);
s+=1;
}
}
/*
for(i=0;i<M;i++)
{
x[i]=X[temp[i]];
fout << "x[" << i << "]=" << x[i] << endl;
}
*/
}
void EPC::transfer2()//把牛顿下山法中的x矩阵转化为牛顿下山法中的a[i][j]矩阵。
{
int s=0;
ofstream fout("transfer2-check.txt");
for(i=1;i<=N;i++)
{
for(j=0;j<=i;j++)
{
a[i][j]=x[s];
s+=1;
}
}
/*
for(i=0;i<=N;i++)
{
for(j=0;j<=N;j++)
fout << "a[" << i << "][" << j << "]=" << a[i][j] << " ";
fout << endl;
}
*/
}
double EPC::factorial(int k)//计算K的二阶乘
{
double sum=1.0;
for(int i=1;i<=k;i++)
{
if((i%2)==0)
sum*=1;
else
sum*=i;
}
return sum;
}
void EPC::Kthmoment()//计算K阶矩
{
ofstream fout("KthMoment.txt");
Sigma1=sqrt(PI*S0/Alpha);
Sigma2=sqrt(PI*S0/Alpha);
for(k=0;k<=3*N;k++)
{
if(k==0)
{
EX1[k]=1;
EX2[k]=1;
}
if((k%2)!=0&&(k!=0))
{
EX1[k]=0;
EX2[k]=0;
}
else
{
EX1[k]=pow(Sigma1,k)*factorial(k-1);
EX2[k]=pow(Sigma2,k)*factorial(k-1);
}
/*
fout << k << "!! = " << factorial(k) << endl << endl;
fout << "EX1[" << k << "]=" << EX1[k] << endl << endl;
fout << "EX2[" << k << "]=" << EX2[k] << endl << endl;
*/
}
}
void EPC::func()
{
ofstream fout("check.txt");
int i,j,c,d,k,l;
int s=0;
for(l=1;l<=N;l++)
{
for(k=0;k<=l;k++)
{
double sum[6]={0,0,0,0,0,0};
for(i=1;i<=N;i++)
{
for(j=0;j<=i;j++)
{
if((i-j-1>=0))
sum[0]+=(a[i][j]*(i-j)*EX1[i-j-1+k]*EX2[j+1+l-k]);
else sum[0]+=0;
if((j-1>=0))
sum[1]+=(Alpha*a[i][j]*j*EX1[i-j+k]*EX2[j+l-k]);
else sum[1]+=0;
if((j-1>=0))
sum[2]+=(a[i][j]*j*EX1[i-j+1+k]*EX2[j-1+l-k]);
else sum[2]+=0;
if((j-2>=0))
sum[3]+=(PI*Alpha*K22*a[i][j]*j*(j-1)*EX1[i-j+k]*EX2[j-2+l-k]);
else sum[3]+=0;
for(c=1;c<=N;c++)
{
for(d=0;d<=c;d++)
{
if((j-1>=0)&&(d-1>=0))
sum[4]+=(a[i][j]*a[c][d]*j*d*EX1[i+c-j-d+k]*EX2[j+d-2+l-k]);
else sum[4]+=0;
}
}
sum[5]+=(Alpha*EX1[k]*EX2[l-k]);
}
}
f[s]=(sum[0]-sum[1]-sum[2]-sum[3]-sum[4]+sum[5]);
s+=1;
}
}
for(s=0;s<M;s++)
{
//cout << "f(" << s << ")=" << f[s] << endl;
fout << "f(" << s << ")=" << f[s] << endl;
}
//cout << "函数值计算完毕" <<endl;
}
void EPC::netn_root () //执行拟牛顿法
{
ofstream fout("check2.txt");
double am,z,beta,d;
l=0; am=1.0+eps;
while (am>=eps)
{
func();
am=0.0;
for (i=0; i<=M-1; i++)
{
B[i] = f[i];
z=fabs(f[i]);
if (z>am) am=z;
}
if (am>=eps)
{
l=l+1;
if (l==Imax)
{
cout <<"\n可能不收敛！" <<endl;
return;
}
transfer1();
for (j=0; j<=M-1; j++)
{
z=x[j]; x[j]=x[j]+h;
func();
for (i=0; i<=M-1; i++) A[i][j] = f[i];
x[j]=z;
}
for(i=0;i<=M-1;i++)
{
for(j=0;j<=M-1;j++)
fout << A[i][j] << " ";
fout << endl;
}
gauss();
beta=1.0;
for (i=0; i<=M-1; i++) beta=beta-B[i];
if (fabs(beta)+1.0==1.0)
{
cout <<"\nBeta=0. 程序工作失败！" <<endl;
return;
}
d=h/beta;
for (i=0; i<=M-1; i++) x[i]=x[i]-d*B[i];
h=t*h;
}
}
cout <<"\n实际迭代次数为: " <<l <<endl;
}
void EPC::gauss () //执行全选主元Gauss消去法
{
int *js,l,k,i,j,is;
double d,t;
js = new int[M];
l=1;
for (k=0; k<=M-2; k++)
{
d=0.0;
for (i=k;i<=M-1;i++)
for (j=k;j<=M-1;j++)
{
t=fabs(A[i][j]);
if (t>d) { d=t; js[k]=j; is=i;}
}
if (d+1.0==1.0) l=0;
else
{ if (js[k]!=k)
for (i=0;i<=M-1;i++)
{
t=A[i][k];
A[i][k]=A[i][js[k]];
A[i][js[k]]=t;
}
if (is!=k)
{
for (j=k;j<=M-1;j++)
{
t=A[k][j];
A[k][j]=A[is][j];
A[is][j]=t;
}
t=B[k]; B[k]=B[is]; B[is]=t;
}
}
if (l==0)
{
delete [] js;
cout <<"\n系数矩阵奇异！无解." <<endl;
return;
}
d=A[k][k];
for (j=k+1;j<=M-1;j++)
A[k][j]=A[k][j]/d;
B[k]=B[k]/d;
for (i=k+1;i<=M-1;i++)
{
for (j=k+1;j<=M-1;j++)
A[i][j]=A[i][j]-A[i][k]*A[k][j];
B[i]=B[i]-A[i][k]*B[k];
}
}
d=A[M-1][M-1];
if (fabs(d)+1.0==1.0)
{
delete [] js;
cout <<"\n系数矩阵奇异！无解." <<endl;
return;
}
B[M-1]=B[M-1]/d;
for (i=M-2;i>=0;i--)
{
t=0.0;
for (j=i+1;j<=M-1;j++)
t=t+A[i][j]*B[j];
B[i]=B[i]-t;
}
js[M-1]=M-1;
for (k=M-1;k>=0;k--)
if (js[k]!=k)
{
t=B[k]; B[k]=B[js[k]]; B[js[k]]=t;
}
delete [] js;
}
void EPC::output () //输出非线性方程组的一组解到文件并显示
{
int k;
ofstream fout ("Results.txt");
if (!fout)
{ cout <<"\n不能打开Result.txt " << endl; exit(1); }
fout <<endl; cout <<endl;
for (k=0; k<M; k++)
{
fout <<x[k] <<" ";
cout <<x[k] <<" ";
}
fout <<endl; cout <<endl;
fout.close ();
}
int _tmain(int argc, _TCHAR* argv[])
{
EPC example1;
example1.input1();
example1.input2();
cout << "基本参数输入完毕" << endl;
//example1.transfer1();
//example1.transfer2();
example1.Kthmoment();
example1.func();
cout << "执行完毕" << endl;
example1.netn_root();
//example1.output();
getch();
return 0;
}
3 个解决方案
#1
因为害怕一次发不完，我在此楼做一些说明
我觉得出问题的是在拟牛顿法中间的各个transfer函数之间的衔接上
因为我的func（）这个函数是用a[i][j]这样一个二维数组作为参数的，但是拟牛顿法又是用x[j]这样一个一维数组作为循环的，所以我分别用了两个transfer函数，一个作用是把二维数组转化为一维数组，另外一个就是把一维数组转化为二维数组。
我觉得出问题的是在拟牛顿法中间的各个transfer函数之间的衔接上
因为我的func（）这个函数是用a[i][j]这样一个二维数组作为参数的，但是拟牛顿法又是用x[j]这样一个一维数组作为循环的，所以我分别用了两个transfer函数，一个作用是把二维数组转化为一维数组，另外一个就是把一维数组转化为二维数组。
#2
学习，。。
#3
早上起来自顶，昨晚想了半夜还是没想出来有可能问题出在哪，求指点求指导啊，前辈高手老师们快出现吧~~！
#1
因为害怕一次发不完，我在此楼做一些说明
我觉得出问题的是在拟牛顿法中间的各个transfer函数之间的衔接上
因为我的func（）这个函数是用a[i][j]这样一个二维数组作为参数的，但是拟牛顿法又是用x[j]这样一个一维数组作为循环的，所以我分别用了两个transfer函数，一个作用是把二维数组转化为一维数组，另外一个就是把一维数组转化为二维数组。
我觉得出问题的是在拟牛顿法中间的各个transfer函数之间的衔接上
因为我的func（）这个函数是用a[i][j]这样一个二维数组作为参数的，但是拟牛顿法又是用x[j]这样一个一维数组作为循环的，所以我分别用了两个transfer函数，一个作用是把二维数组转化为一维数组，另外一个就是把一维数组转化为二维数组。
#2
学习，。。
#3
早上起来自顶，昨晚想了半夜还是没想出来有可能问题出在哪，求指点求指导啊，前辈高手老师们快出现吧~~！