?? nlequation.cpp
字號:
}
}
}
*x=xx;
}
//////////////////////////////////////////////////////////////////////
// 內部函數(shù)
//////////////////////////////////////////////////////////////////////
double CNLequation::rnd(double* r)
{
int m;
double s,u,v,p;
s=65536.0;
u=2053.0;
v=13849.0;
m=(int)(*r/s);
*r=*r-m*s;
*r=u*(*r)+v;
m=(int)(*r/s);
*r=*r-m*s; p=*r/s;
return(p);
}
//////////////////////////////////////////////////////////////////////
// 求實函數(shù)或復函數(shù)方程一個復根的蒙特卡洛法
//
// 調用時，須覆蓋計算方程左端函數(shù)的模值||f(x, y)||的虛函數(shù)
// double Func(double x, double y)
//
// 參數(shù)：
// 1. double* x - 傳入初值（猜測解）的實部，返回求得的根的實部
// 2. double* y - 傳入初值（猜測解）的虛部，返回求得的根的虛部
// 3. double xStart - 均勻分布的端點初值
// 4. int nControlB - 控制參數(shù)
// 5. double eps - 控制精度，默認值為0.000001
// 6. double xi[] - 一維數(shù)組，長度為n，返回n個根的虛部
//
// 返回值：無
//////////////////////////////////////////////////////////////////////
void CNLequation::GetRootMonteCarlo(double* x, double* y, double xStart, int nControlB, double eps /*= 0.000001*/)
{
int k;
double xx,yy,a,r,z,x1,y1,z1;
// 求解條件與初值
a=xStart;
k=1;
r=1.0;
xx=*x;
yy=*y;
z=Func(xx,yy);
// 精度控制求解
while (a>=eps)
{
x1=-a+2.0*a*rnd(&r);
x1=xx+x1;
y1=-a+2.0*a*rnd(&r);
y1=yy+y1;
z1=Func(x1,y1);
k=k+1;
if (z1>=z)
{
if (k>nControlB)
{
k=1;
a=a/2.0;
}
}
else
{
k=1;
xx=x1;
yy=y1;
z=z1;
if (z<eps)
{
*x=xx;
*y=yy;
return;
}
}
}
*x=xx;
*y=yy;
}
//////////////////////////////////////////////////////////////////////
// 求非線性方程組一組實根的梯度法
//
// 調用時，須覆蓋計算方程左端函數(shù)f(x)值及其偏導數(shù)值的虛函數(shù)
// double Func(double x[], double y[])
//
// 參數(shù)：
// 1. int n - 方程的個數(shù)，也是未知數(shù)的個數(shù)
// 2. double x[] - 一維數(shù)組，長度為n，存放一組初值x0, x1, …, xn-1，
// 返回時存放方程組的一組實根
// 3. int nMaxIt - 迭代次數(shù)，默認值為60
// 4. double eps - 控制精度，默認值為0.000001
//
// 返回值：BOOL 型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CNLequation::GetRootsetGrad(int n, double x[], int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
{
int l,j;
double f,d,s,*y;
y = new double[n];
l=nMaxIt;
f=Func(x,y);
// 控制精度，迭代求解
while (f>=eps)
{
l=l-1;
if (l==0)
{
delete[] y;
return TRUE;
}
d=0.0;
for (j=0; j<=n-1; j++)
d=d+y[j]*y[j];
if (d+1.0==1.0)
{
delete[] y;
return FALSE;
}
s=f/d;
for (j=0; j<=n-1; j++)
x[j]=x[j]-s*y[j];
f=Func(x,y);
}
delete[] y;
// 是否在有效迭代次數(shù)內達到精度
return (nMaxIt>l);
}
//////////////////////////////////////////////////////////////////////
// 求非線性方程組一組實根的擬牛頓法
//
// 調用時，須覆蓋計算方程左端函數(shù)f(x)值及其偏導數(shù)值的虛函數(shù)
// double Func(double x[], double y[])
//
// 參數(shù)：
// 1. int n - 方程的個數(shù)，也是未知數(shù)的個數(shù)
// 2. double x[] - 一維數(shù)組，長度為n，存放一組初值x0, x1, …, xn-1，
// 返回時存放方程組的一組實根
// 3. double t - 控制h大小的變量，0<t<1
// 4. double h - 增量初值
// 3. int nMaxIt - 迭代次數(shù)，默認值為60
// 4. double eps - 控制精度，默認值為0.000001
//
// 返回值：BOOL 型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CNLequation::GetRootsetNewton(int n, double x[], double t, double h, int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
{
int i,j,l;
double am,z,beta,d,*y,*a,*b;
y = new double[n];
// 構造矩陣
CMatrix mtxCoef(n, n);
CMatrix mtxConst(n, 1);
a = mtxCoef.GetData();
b = mtxConst.GetData();
// 迭代求解
l=nMaxIt;
am=1.0+eps;
while (am>=eps)
{
Func(x,b);
am=0.0;
for (i=0; i<=n-1; i++)
{
z=fabs(b[i]);
if (z>am)
am=z;
}
if (am>=eps)
{
l=l-1;
if (l==0)
{
delete[] y;
return FALSE;
}
for (j=0; j<=n-1; j++)
{
z=x[j];
x[j]=x[j]+h;
Func(x,y);
for (i=0; i<=n-1; i++)
a[i*n+j]=y[i];
x[j]=z;
}
// 調用全選主元高斯消元法
CLEquations leqs(mtxCoef, mtxConst);
CMatrix mtxResult;
if (! leqs.GetRootsetGauss(mtxResult))
{
delete[] y;
return FALSE;
}
mtxConst = mtxResult;
b = mtxConst.GetData();
beta=1.0;
for (i=0; i<=n-1; i++)
beta=beta-b[i];
if (beta == 0.0)
{
delete[] y;
return FALSE;
}
d=h/beta;
for (i=0; i<=n-1; i++)
x[i]=x[i]-d*b[i];
h=t*h;
}
}
delete[] y;
// 是否在有效迭代次數(shù)內達到精度
return (nMaxIt>l);
}
//////////////////////////////////////////////////////////////////////
// 求非線性方程組最小二乘解的廣義逆法
//
// 調用時，1. 須覆蓋計算方程左端函數(shù)f(x)值及其偏導數(shù)值的虛函數(shù)
// double Func(double x[], double y[])
// 2. 須覆蓋計算雅可比矩陣函數(shù)的虛函數(shù)
// double FuncMJ(int n, double x[], double y[])
//
// 參數(shù)：
// 1. int m - 方程的個數(shù)
// 2. int n - 未知數(shù)的個數(shù)
// 3. double x[] - 一維數(shù)組，長度為n，存放一組初值x0, x1, …, xn-1，要求
// 不全為0，返回時存放方程組的最小二乘解，當m=n時，即是
// 非線性方程組的解
// 4. double eps1 - 最小二乘解的精度控制精度，默認值為0.000001
// 5. double eps2 - 奇異值分解的精度控制精度，默認值為0.000001
//
// 返回值：BOOL 型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CNLequation::GetRootsetGinv(int m, int n, double x[], double eps1 /*= 0.000001*/, double eps2 /*= 0.000001*/)
{
int i,j,k,l,kk,jt;
double y[10],b[10],alpha,z,h2,y1,y2,y3,y0,h1;
double *p,*d,*dx,*w;
// 控制參數(shù)
int ka = max(m, n)+1;
w = new double[ka];
// 設定迭代次數(shù)為60，迭代求解
l=60;
alpha=1.0;
while (l>0)
{
CMatrix mtxP(m, n), mtxD(m, 1);
p = mtxP.GetData();
d = mtxD.GetData();
Func(x,d);
FuncMJ(n,x,p);
// 構造線性方程組
CLEquations leqs(mtxP, mtxD);
// 臨時矩陣
CMatrix mtxAP, mtxU, mtxV;
// 解矩陣
CMatrix mtxDX;
// 基于廣義逆的最小二乘解
if (! leqs.GetRootsetGinv(mtxDX, mtxAP, mtxU, mtxV, eps2))
{
delete[] w;
return FALSE;
}
dx = mtxDX.GetData();
j=0;
jt=1;
h2=0.0;
while (jt==1)
{
jt=0;
if (j<=2)
z=alpha+0.01*j;
else
z=h2;
for (i=0; i<=n-1; i++)
w[i]=x[i]-z*dx[i];
Func(w,d);
y1=0.0;
for (i=0; i<=m-1; i++)
y1=y1+d[i]*d[i];
for (i=0; i<=n-1; i++)
w[i]=x[i]-(z+0.00001)*dx[i];
Func(w,d);
y2=0.0;
for (i=0; i<=m-1; i++)
y2=y2+d[i]*d[i];
y0=(y2-y1)/0.00001;
if (fabs(y0)>1.0e-10)
{
h1=y0; h2=z;
if (j==0)
{
y[0]=h1;
b[0]=h2;
}
else
{
y[j]=h1;
kk=0;
k=0;
while ((kk==0)&&(k<=j-1))
{
y3=h2-b[k];
if (fabs(y3)+1.0==1.0)
kk=1;
else
h2=(h1-y[k])/y3;
k=k+1;
}
b[j]=h2;
if (kk!=0)
b[j]=1.0e+35;
h2=0.0;
for (k=j-1; k>=0; k--)
h2=-y[k]/(b[k+1]+h2);
h2=h2+b[0];
}
j=j+1;
if (j<=7)
jt=1;
else
z=h2;
}
}
alpha=z;
y1=0.0;
y2=0.0;
for (i=0; i<=n-1; i++)
{
dx[i]=-alpha*dx[i];
x[i]=x[i]+dx[i];
y1=y1+fabs(dx[i]);
y2=y2+fabs(x[i]);
}
// 求解成功
if (y1<eps1*y2)
{
delete[] w;
return TRUE;
}
l=l-1;
}
delete[] w;
// 求解失敗
return FALSE;
}
//////////////////////////////////////////////////////////////////////
// 求非線性方程組一組實根的蒙特卡洛法
//
// 調用時，須覆蓋計算方程左端模函數(shù)值||F||的虛函數(shù)
// double Func(int n, double x[])
// 其返回值為Sqr(f1*f1 + f2*f2 + … + fn*fn)
//
// 參數(shù)：
// 1. double x[] - 一維數(shù)組，長度為n，存放一組初值，返回時存放方程的一組實根
// 2. double xStart - 均勻分布的端點初值
// 3. int nControlB - 控制參數(shù)
// 4. double eps - 控制精度，默認值為0.000001
//
// 返回值：無
//////////////////////////////////////////////////////////////////////
void CNLequation::GetRootsetMonteCarlo(int n, double x[], double xStart, int nControlB, double eps /*= 0.000001*/)
{
int k,i;
double a,r,*y,z,z1;
y = new double[n];
// 初值
a=xStart;
k=1;
r=1.0;
z = Func(n, x);
// 用精度控制迭代求解
while (a>=eps)
{
for (i=0; i<=n-1; i++)
y[i]=-a+2.0*a*rnd(&r)+x[i];
z1=Func(n,y);
k=k+1;
if (z1>=z)
{
if (k>nControlB)
{
k=1;
a=a/2.0;
}
}
else
{
k=1;
for (i=0; i<=n-1; i++)
x[i]=y[i];
// 求解成功
z=z1;
if (z<eps)
{
delete[] y;
return;
}
}
}
delete[] y;
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -