?? nlequation.cpp
字號:
//////////////////////////////////////////////////////////////////////
// NLequation.cpp
//
// 求解線性方程組的類 CNLequation 的實現代碼
//
// 周長發編制, 2002/8
//////////////////////////////////////////////////////////////////////
#include "stdafx.h"
#include "NLequation.h"
#include "LEquations.h"
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////
// 基本構造函數
//////////////////////////////////////////////////////////////////////
CNLequation::CNLequation()
{
}
//////////////////////////////////////////////////////////////////////
// 析構函數
//////////////////////////////////////////////////////////////////////
CNLequation::~CNLequation()
{
}
//////////////////////////////////////////////////////////////////////
// 求非線性方程實根的對分法
//
// 調用時，須覆蓋計算方程左端函數f(x)值的虛函數double Func(double x)
//
// 參數：
// 1. int nNumRoots - 在[xStart, xEnd]內實根個數的預估值
// 2. double x[] - 一維數組，長度為m。返回在區間[xStart, xEnd]內搜索到
// 的實根，實根個數由函數值返回
// 3. double xStart - 求根區間的左端點
// 4. double xEnd - 求根區間的右端點
// 5. double dblStep - 搜索求根時采用的步長
// 6. double eps - 精度控制參數，默認值為0.000001
//
// 返回值：int 型，求得的實根的數目
//////////////////////////////////////////////////////////////////////
int CNLequation::GetRootBisect(int nNumRoots, double x[], double xStart, double xEnd, double dblStep, double eps /*= 0.000001*/)
{
int n,js;
double z,y,z1,y1,z0,y0;
// 根的個數清0
n = 0;
// 從左端點開始搜索
z = xStart;
y = Func(z);
// 循環求解
while ((z<=xEnd+dblStep/2.0)&&(n!=nNumRoots))
{
if (fabs(y)<eps)
{
n=n+1;
x[n-1]=z;
z=z+dblStep/2.0;
y=Func(z);
}
else
{
z1=z+dblStep;
y1=Func(z1);
if (fabs(y1)<eps)
{
n=n+1;
x[n-1]=z1;
z=z1+dblStep/2.0;
y=Func(z);
}
else if (y*y1>0.0)
{
y=y1;
z=z1;
}
else
{
js=0;
while (js==0)
{
if (fabs(z1-z)<eps)
{
n=n+1;
x[n-1]=(z1+z)/2.0;
z=z1+dblStep/2.0; y=Func(z);
js=1;
}
else
{
z0=(z1+z)/2.0;
y0=Func(z0);
if (fabs(y0)<eps)
{
x[n]=z0;
n=n+1;
js=1;
z=z0+dblStep/2.0;
y=Func(z);
}
else if ((y*y0)<0.0)
{
z1=z0;
y1=y0;
}
else
{
z=z0;
y=y0;
}
}
}
}
}
}
// 返回實根的數目
return(n);
}
//////////////////////////////////////////////////////////////////////
// 求非線性方程一個實根的牛頓法
//
// 調用時，須覆蓋計算方程左端函數f(x)及其一階導數f'(x)值的虛函數
// void Func(double x, double y[])
// y(0) 返回f(x)的值
// y(1) 返回f'(x)的值
//
// 參數：
// 1. double *x - 傳入迭代初值（猜測解），返回在區間求得的一個實根
// 2. int nMaxIt - 遞歸次數，默認值為60
// 3. double eps - 精度控制參數，默認值為0.000001
//
// 返回值：BOOL 型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CNLequation::GetRootNewton(double* x, int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
{
int l;
double y[2],d,p,x0,x1;
// 條件值
l=nMaxIt;
x0=*x;
Func(x0,y);
// 求解，控制精度
d=eps+1.0;
while ((d>=eps)&&(l!=0))
{
if (y[1] == 0.0)
return FALSE;
x1=x0-y[0]/y[1];
Func(x1,y);
d=fabs(x1-x0);
p=fabs(y[0]);
if (p>d)
d=p;
x0=x1;
l=l-1;
}
*x=x1;
return TRUE;
}
//////////////////////////////////////////////////////////////////////
// 求非線性方程一個實根的埃特金迭代法
//
// 調用時，須覆蓋計算方程左端函數f(x)值的虛函數double Func(double x)
//
// 參數：
// 1. double *x - 傳入迭代初值（猜測解），返回在區間求得的一個實根
// 2. int nMaxIt - 遞歸次數，默認值為60
// 3. double eps - 精度控制參數，默認值為0.000001
//
// 返回值：BOOL 型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CNLequation::GetRootAitken(double* x, int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
{
int flag,l;
double u,v,x0;
// 求解條件
l=0;
x0=*x;
flag=0;
// 迭代求解
while ((flag==0)&&(l!=nMaxIt))
{
l=l+1;
u=Func(x0);
v=Func(u);
if (fabs(u-v)<eps)
{
x0=v;
flag=1;
}
else
x0=v-(v-u)*(v-u)/(v-2.0*u+x0);
}
*x=x0;
// 是否在指定的迭代次數內達到求解精度
return (nMaxIt > l);
}
//////////////////////////////////////////////////////////////////////
// 求非線性方程一個實根的連分式解法
//
// 調用時，須覆蓋計算方程左端函數f(x)值的虛函數double Func(double x)
//
// 參數：
// 1. double *x - 傳入迭代初值（猜測解），返回在區間求得的一個實根
// 2. double eps - 精度控制參數，默認值為0.000001
//
// 返回值：BOOL 型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CNLequation::GetRootPq(double* x, double eps /*= 0.000001*/)
{
int i,j,m,it,l;
double a[10],y[10],z,h,x0,q;
// 求解條件
l=10;
q=1.0e+35;
x0=*x;
h=0.0;
// 連分式求解
while (l!=0)
{
l=l-1;
j=0;
it=l;
while (j<=7)
{
if (j<=2)
z=x0+0.1*j;
else
z=h;
y[j]=Func(z);
h=z;
if (j==0)
a[0]=z;
else
{
m=0;
i=0;
while ((m==0)&&(i<=j-1))
{
if (fabs(h-a[i])+1.0==1.0)
m=1;
else
h=(y[j]-y[i])/(h-a[i]);
i=i+1;
}
a[j]=h;
if (m!=0)
a[j]=q;
h=0.0;
for (i=j-1; i>=0; i--)
{
if (fabs(a[i+1]+h)+1.0==1.0)
h=q;
else
h=-y[i]/(a[i+1]+h);
}
h=h+a[0];
}
if (fabs(y[j])>=eps)
j=j+1;
else
{
j=10;
l=0;
}
}
x0=h;
}
*x=h;
// 是否在10階連分式內求的實根？
return (10>it);
}
//////////////////////////////////////////////////////////////////////
// 求實系數代數方程全部根的QR方法
//
// 參數：
// 1. int n - 多項式方程的次數
// 2. double dblCoef[] - 一維數組，長度為n+1，按降冪次序依次存放n次多項式方程的n+1個系數
// 3. double xr[] - 一維數組，長度為n，返回n個根的實部
// 4. double xi[] - 一維數組，長度為n，返回n個根的虛部
// 5. int nMaxIt - 迭代次數，默認值為 60
// 6. double eps - 精度控制參數，默認值為 0.000001
//
// 返回值：BOOL 型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CNLequation::GetRootQr(int n, double dblCoef[], double xr[], double xi[], int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
{
// 初始化矩陣
CMatrix mtxQ;
mtxQ.Init(n, n);
double *q = mtxQ.GetData();
// 構造赫申伯格矩陣
for (int j=0; j<=n-1; j++)
q[j]=-dblCoef[n-j-1]/dblCoef[n];
for (j=n; j<=n*n-1; j++)
q[j]=0.0;
for (int i=0; i<=n-2; i++)
q[(i+1)*n+i]=1.0;
// 求赫申伯格矩陣的特征值和特征向量，即為方程的解
if (mtxQ.HBergEigenv(xr, xi, nMaxIt, eps))
return TRUE;
return FALSE;
}
//////////////////////////////////////////////////////////////////////
// 求實系數代數方程全部根的牛頓下山法
//
// 參數：
// 1. int n - 多項式方程的次數
// 2. double dblCoef[] - 一維數組，長度為n+1，按降冪次序依次存放n次多項式方程的n+1個系數
// 3. double xr[] - 一維數組，長度為n，返回n個根的實部
// 4. double xi[] - 一維數組，長度為n，返回n個根的虛部
//
// 返回值：BOOL 型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CNLequation::GetRootNewtonDownHill(int n, double dblCoef[], double xr[], double xi[])
{
int m,i,jt,k,is,it;
double t,x,y,x1,y1,dx,dy,p,q,w,dd,dc,c;
double g,u,v,pq,g1,u1,v1;
// 初始判斷
m=n;
while ((m>0)&&(fabs(dblCoef[m])+1.0==1.0))
m=m-1;
// 求解失敗
if (m<=0)
return FALSE;
for (i=0; i<=m; i++)
dblCoef[i]=dblCoef[i]/dblCoef[m];
for (i=0; i<=m/2; i++)
{
w=dblCoef[i];
dblCoef[i]=dblCoef[m-i];
dblCoef[m-i]=w;
}
// 迭代求解
k=m;
is=0;
w=1.0;
jt=1;
while (jt==1)
{
pq=fabs(dblCoef[k]);
while (pq<1.0e-12)
{
xr[k-1]=0.0;
xi[k-1]=0.0;
k=k-1;
if (k==1)
{
xr[0]=-dblCoef[1]*w/dblCoef[0];
xi[0]=0.0;
return TRUE;
}
pq=fabs(dblCoef[k]);
}
q=log(pq);
q=q/(1.0*k);
q=exp(q);
p=q;
w=w*p;
for (i=1; i<=k; i++)
{
dblCoef[i]=dblCoef[i]/q;
q=q*p;
}
x=0.0001;
x1=x;
y=0.2;
y1=y;
dx=1.0;
g=1.0e+37;
l40:
u=dblCoef[0]; v=0.0;
for (i=1; i<=k; i++)
{
p=u*x1;
q=v*y1;
pq=(u+v)*(x1+y1);
u=p-q+dblCoef[i];
v=pq-p-q;
}
g1=u*u+v*v;
if (g1>=g)
{
if (is!=0)
{
it=1;
g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
if (it==0)
goto l40;
}
else
{
g60(&t,&x,&y,&x1,&y1,&dx,&dy,&p,&q,&k,&it);
if (t>=1.0e-03)
goto l40;
if (g>1.0e-18)
{
it=0;
g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
if (it==0)
goto l40;
}
}
g90(xr,xi,dblCoef,&x,&y,&p,&q,&w,&k);
}
else
{
g=g1;
x=x1;
y=y1;
is=0;
if (g<=1.0e-22)
g90(xr,xi,dblCoef,&x,&y,&p,&q,&w,&k);
else
{
u1=k*dblCoef[0];
v1=0.0;
for (i=2; i<=k; i++)
{
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -