?? nlequations.cs
字號:
/*
* 求解非線性方程組的類 NLEquations
*
* 周長發編制
*/
using System;
namespace CSharpAlgorithm.Algorithm
{
/**
* 求解非線性方程組的類 NLEquations
*
* @author 周長發
* @version 1.0
*/
public abstract class NLEquations
{
/**
* 虛函數：計算方程左端函數值，必須在引申類中覆蓋該類函數
*
* @param x - 變量
* @return 函數值
*/
protected virtual double Func(double x)
{
return 0;
}
/**
* 虛函數：計算方程左端函數值，必須在引申類中覆蓋該類函數
*
* @param x - 變量值數組
* @return 函數值
*/
protected virtual double Func(double[] x)
{
return 0;
}
/**
* 虛函數：計算方程左端函數值，必須在引申類中覆蓋該類函數
*
* @param x - 變量
* @param y - 函數值數組
*/
protected virtual void Func(double x, double[] y)
{
}
/**
* 虛函數：計算方程左端函數值，必須在引申類中覆蓋該類函數
*
* @param x - 二元函數的變量
* @param y - 二元函數的變量
* @return 函數值
*/
protected virtual double Func(double x, double y)
{
return 0;
}
/**
* 虛函數：計算方程左端函數值，必須在引申類中覆蓋該類函數
*
* @param x - 二元函數的變量值數組
* @param y - 二元函數的變量值數組
* @return 函數值
*/
protected virtual double Func(double[] x, double[] y)
{
return 0;
}
/**
* 虛函數：計算方程左端函數值，必須在引申類中覆蓋該類函數
*
* @param x - 已知變量值數組
* @param p - 已知函數值數組
*/
protected virtual void FuncMJ(double[] x, double[] p)
{
}
/**
* 基本構造函數
*/
public NLEquations()
{
}
/**
* 求非線性方程實根的對分法
*
* 調用時，須覆蓋計算方程左端函數f(x)值的虛函數
* double Func(double x)
*
* @param nNumRoots - 在[xStart, xEnd]內實根個數的預估值
* @param x - 一維數組，長度為m。返回在區間[xStart, xEnd]內搜索到的實根，
* 實根個數由函數值返回
* @param xStart - 求根區間的左端點
* @param xEnd - 求根區間的右端點
* @param dblStep - 搜索求根時采用的步長
* @param eps - 精度控制參數
* @return int 型，求得的實根的數目
*/
public int GetRootBisect(int nNumRoots, double[] x, double xStart, double xEnd, double dblStep, double eps)
{
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 (Math.Abs(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 (Math.Abs(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 (Math.Abs(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 (Math.Abs(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)的值
*
* @param x - 傳入迭代初值（猜測解），返回在區間求得的一個實根
* @param nMaxIt - 遞歸次數
* @param eps - 精度控制參數
* @return bool 型，求解是否成功
*/
public bool GetRootNewton(ref double x, int nMaxIt, double eps)
{
int l;
double d,p,x0,x1=0.0;
double[] y = new double[2];
// 條件值
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=Math.Abs(x1-x0);
p=Math.Abs(y[0]);
if (p>d)
d=p;
x0=x1;
l=l-1;
}
x = x1;
return true;
}
/**
* 求非線性方程一個實根的埃特金迭代法
*
* 調用時，須覆蓋計算方程左端函數f(x)值的虛函數
* double Func(double x)
*
* @param x - 傳入迭代初值（猜測解），返回在區間求得的一個實根
* @param nMaxIt - 遞歸次數
* @param eps - 精度控制參數
* @return bool 型，求解是否成功
*/
public bool GetRootAitken(ref double x, int nMaxIt, double eps)
{
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 (Math.Abs(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)
*
* @param x - 傳入迭代初值（猜測解），返回在區間求得的一個實根
* @param eps - 精度控制參數
* @return bool 型，求解是否成功
*/
public bool GetRootPq(ref double x, double eps)
{
int i,j,m,it=0,l;
double z,h,x0,q;
double[] a = new double[10];
double[] y = new double[10];
// 求解條件
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 (Math.Abs(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 (Math.Abs(a[i+1]+h)+1.0==1.0)
h=q;
else
h=-y[i]/(a[i+1]+h);
}
h=h+a[0];
}
if (Math.Abs(y[j])>=eps)
j=j+1;
else
{
j=10;
l=0;
}
}
x0=h;
}
x = h;
// 是否在10階連分式內求的實根？
return (10>it);
}
/**
* 求實系數代數方程全部根的QR方法
*
* @param n - 多項式方程的次數
* @param dblCoef - 一維數組，長度為n+1，按降冪次序依次存放n次多項式方程的
* n+1個系數
* @param xr - 一維數組，長度為n，返回n個根的實部
* @param xi - 一維數組，長度為n，返回n個根的虛部
* @param nMaxIt - 迭代次數
* @param eps - 精度控制參數
* @return bool 型，求解是否成功
*/
public bool GetRootQr(int n, double[] dblCoef, double[] xr, double[] xi, int nMaxIt, double eps)
{
// 初始化矩陣
Matrix mtxQ = new Matrix();
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 (int 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.ComputeEvHBerg(xr, xi, nMaxIt, eps))
return true;
return false;
}
/**
* 求實系數代數方程全部根的牛頓下山法
*
* @param n - 多項式方程的次數
* @param dblCoef - 一維數組，長度為n+1，按降冪次序依次存放n次多項式方程的
* n+1個系數
* @param xr - 一維數組，長度為n，返回n個根的實部
* @param xi - 一維數組，長度為n，返回n個根的虛部
* @return bool 型，求解是否成功
*/
public bool GetRootNewtonDownHill(int n, double[] dblCoef, double[] xr, double[] xi)
{
int m=0,i=0,jt=0,k=0,nis=0,it=0;
double t=0,x=0,y=0,x1=0,y1=0,dx=0,dy=0,p=0,q=0,w=0,dd=0,dc=0,c=0;
double g=0,u=0,v=0,pq=0,g1=0,u1=0,v1=0;
// 初始判斷
m=n;
while ((m>0)&&(Math.Abs(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;
nis=0;
w=1.0;
jt=1;
while (jt==1)
{
pq=Math.Abs(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=Math.Abs(dblCoef[k]);
}
q=Math.Log(pq);
q=q/(1.0*k);
q=Math.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;
while (true)
{
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 (nis!=0)
{
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -