?? nonlinearequation.inl
字號:
//NonLinearEquation.inl 非線性方程(組)求解函數(shù)(方法)定義
// Ver 1.0.0.0
// 版權(quán)所有(C) 何渝, 2002
// 最后修改: 2002.5.31.
#ifndef _NONLINEAREQUATION_INL
#define _NONLINEAREQUATION_INL
//用二分法搜索方程f(x)=0在區(qū)間[a,b]內(nèi)的全部實根
template <class _Ty>
inline size_t
RootHalves(_Ty a, _Ty b, _Ty step, _Ty eps, valarray<_Ty>& x, size_t m)
{
int n(0), js;
_Ty z=a, z1, y, y1, z0, y0;
y = FunctionValueRH(z);
while((z<=b+step/2.0)&&(n!=m))
{
if(Abs(y) < eps)
{
n = n + 1;
x[n-1] = z;
z = z + step / 2.0;
y=FunctionValueRH(z);
}
else
{
z1=z+step;
y1=FunctionValueRH(z1);
if(Abs(y1)<eps)
{
n=n+1;
x[n-1]=z1;
z=z1+step/2.0;
y=FunctionValueRH(z);
}
else
if(y*y1>0.0)
{
y=y1;
z=z1;
}
else
{
js=0;
while(js==0)
{
if(Abs(z1-z)<eps)
{
n=n+1;
x[n-1]=(z1+z)/2.0;
z=z1+step/2.0;
y=FunctionValueRH(z);
js=1;
}
else
{
z0=(z1+z)/2.0;
y0=FunctionValueRH(z0);
if(Abs(y0)<eps)
{
x[n]=z0; n=n+1; js=1;
z=z0+step/2.0;
y=FunctionValueRH(z);
}
else
if((y*y0) < 0.0)
{
z1 = z0;
y1 = y0;
}
else
{
z = z0;
y = y0;
}
}
}
}
}
}
return(n); //返回搜索到的實根個數(shù)
}
//牛頓(Newton)法求解非線性方程一個實根
template <class _Ty>
inline int
RootNewton( _Ty& x, _Ty eps, size_t js)
{
int k = 1, r = js;
_Ty x1;
_Ty x0 = x;
valarray<_Ty> y(2);
FunctionValueRN(x0, y);
_Ty d=eps + 1.0;
while((d > eps || FloatEqual(d, eps))&&(r != 0))
{
if(FloatEqual(y[1],0)) return(-1);
x1 = x0 - y[0] / y[1];
FunctionValueRN(x1,y);
d=Abs(x1 - x0);
_Ty p=Abs(y[0]);
if(p > d) d = p;
x0 = x1;
r -= 1;
}
x = x1;
k = js - r;
return(k);
}
//埃特金(Aitken)法求解非線性方程一個實根
template <class _Ty>
inline int
RootAitken( _Ty& x, _Ty eps, size_t js)
{
_Ty u, v, x0;
int r(0), flag(0);
x0 = x;
while((flag == 0) && (r != js))
{
r++;
u = FunctionValueRA(x0);
v = FunctionValueRA(u);
if(Abs(u-v)<eps)
{
x0 = v;
flag = 1;
}
else x0=v-(v-u)*(v-u)/(v-2.0*u+x0);
}
x = x0;
r = js - r;
return(r);
}
//連分式(Fraction)法求解非線性方程一個實根
template <class _Ty>
inline int
RootFraction(_Ty& x, _Ty eps)
{
int r(10), i, j, m, it;
_Ty x0, q(1.0e+35), h(0), z;
valarray<_Ty> a(10), y(10);
x0 = x;
while(r!=0)
{
r = r - 1;
j = 0;
it = r;
while(j<=7)
{
if(j<=2)
z = x0 + 0.1 * j;
else
z = h;
y[j] = FunctionValueRF(z);
h = z;
if(j==0)
a[0] = z;
else
{
m = 0;
i = 0;
while((m==0)&&(i<j))
{
if(FloatEqual((h-a[i]),0))
m = 1;
else
h=(y[j]-y[i])/(h-a[i]);
i++;
}
a[j] = h;
if(m!=0) a[j]=q;
h=0.0;
for(i=j-1; i>=0; i--)
{
if(FloatEqual((a[i+1]+h),0))
h = q;
else
h=-y[i]/(a[i+1]+h);
}
h = h + a[0];
}
if(Abs(y[j])>=eps) j++;
else
{
j = 10;
r = 0;
}
}
x0 = h;
}
x = h;
return(10-it);
}
//QR法求代數(shù)方程全部根
template <class _Ty, class _Tz>
inline int
RootQR(valarray<_Ty>& a, valarray<_Tz>& x, _Ty eps, size_t jt)
{
int i, j, n;
n = x.size(); //多項式方程的次數(shù)
matrix<_Ty> q(n, n);
for(j=0; j<n; j++)
q(0,j) = -a[n-j-1] / a[n];
for(j=1; j<n; j++)
for(i=0; i<n; i++)
q(j,i) = 0;
for(i=0; i<n-1; i++)
q(i+1, i) = 1;
i = EigenvalueVectorHessenbergQR(q,x,eps,jt); //求全部特征值的QR方法
return(i);
}
//牛頓下山(NewtonHillDown)法求解實系數(shù)代數(shù)方程全部根(實根和復(fù)根)
template <class _Ty, class _Tz>
inline int
RootNewtonHillDown(valarray<_Ty>& a, valarray<_Tz>& cx)
{
int m, i, jt, k, is, it, n;
_Ty t, x, y, x1, y1, dx, dy, dd, dc, c, g, g1;
_Ty w, u, p, q, pq, v, u1, v1;
n = cx.size(); //多項式方程的次數(shù)
m = n;
while((m > 0) && (FloatEqual(a[m], 0.0))) m--;
if(m <= 0) return(-1);
for(i = 0; i <= m; i++) a[i] /= a[m];
for(i = 0; i <= m / 2; i++)
{
w = a[i];
a[i] = a[m - i];
a[m - i] = w;
}
k = m;
is = 0;
w = 1.0;
jt = 1;
while(jt == 1)
{
pq = Abs(a[k]);
while(pq < LONGDOUBLEERROR)
{
cx[k - 1] = complex<_Ty>(0.0, 0.0);
k = k - 1;
if(k == 1)
{
cx[0] = complex<_Ty>(-a[1] * w / a[0], 0.0);
return(1);
}
pq = Abs(a[k]);
}
q = log(pq);
q = q / (1.0 * k);
q = exp(q);
p = q;
w = w * p;
for(i = 1; i <= k; i++)
{
a[i] /= q;
q *= p;
}
x = 0.0001;
x1 = x;
y = 0.2;
y1 = y;
dx = 1.0;
g = 1.0e+37;
zjn: u = a[0];
v = 0.0;
for(i = 1; i <= k; i++)
{
p = u * x1;
q = v * y1;
pq = (u + v) * (x1 + y1);
u = p - q + a[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 zjn;
}
else
{
g60(t, x, y, x1, y1, dx, dy, p, q, k, it);
if(t >= 1.0e-03) goto zjn;
if(g > 1.0e-18)
{
it = 0;
g65(x, y, x1, y1, dx, dy, dd, dc, c, k, is, it);
if(it == 0) goto zjn;
}
}
g90(cx , a, x, y, p, q, w, k);
}
else
{
g = g1;
x = x1;
y = y1;
is = 0;
if(g <= 1.0e-22)
g90(cx, a, x, y, p, q, w, k);
else
{
u1 = k * a[0];
v1 = 0.0;
for(i = 2; i <= k; i++)
{
p = u1 * x;
q = v1 * y;
pq = (u1 + v1) * (x + y);
u1 = p - q + (k - i + 1) * a[i - 1];
v1 = pq - p - q;
}
p = u1 * u1 + v1 * v1;
if(p <= 1.0e-20)
{
it = 0;
g65(x, y, x1, y1, dx, dy, dd, dc, c, k, is, it);
if(it == 0) goto zjn;
g90(cx, a, x, y, p, q, w, k);
}
else
{
dx = (u * u1 + v * v1) / p;
dy = (u1 * v - v1 * u) / p;
t = 1.0 + 4.0 / k;
g60(t, x, y, x1, y1, dx, dy, p, q, k, it);
if(t >= 1.0e-03) goto zjn;
if(g > 1.0e-18)
{
it = 0;
g65(x, y, x1, y1, dx, dy, dd, dc, c, k, is, it);
if(it == 0) goto zjn;
}
g90(cx, a, x, y, p, q, w, k);
}
}
}
if(k == 1) jt = 0;
else jt = 1;
}
return(1);
}
//牛頓下山(NewtonHillDown)法求解復(fù)系數(shù)代數(shù)方程全部根(實根和復(fù)根)
//重載RootNewtonHillDown()
template <class _Ty>
inline int
RootNewtonHillDown(complex<_Ty> a[], valarray<complex<_Ty> >& cx)
{
int m, i, jt, k, is, it, n;
_Ty t, x, y, x1, y1, dx, dy, dd, dc, c, g, g1, mo, sb, xb;
_Ty w, u, p, q, pq, v, u1, v1;
n = cx.size(); //多項式方程的次數(shù)
m = n;
while((m > 0) && (FloatEqual(Abs(a[m]), 0))) m--;
if(m <= 0) return(-1);
mo = Abs(a[m]);
for(i = 0; i <= m; i++) a[i] /= mo;
for(i = 0; i <= m / 2; i++)
{
swap(a[i], a[m-i]);
}
k = m;
is = 0;
w = 1.0;
jt = 1;
while(jt == 1)
{
pq = Abs(a[k]);
while(pq < LONGDOUBLEERROR)
{
cx[k - 1] = complex<_Ty>(0.0, 0.0);
k--;
if(k == 1)
{
mo = a[0].real() * a[0].real() + a[0].imag() * a[0].imag();
sb = -w * (a[0].real() * a[1].real() + a[0].imag() * a[1].imag()) / mo;
xb = w * (a[1].real() * a[0].real() - a[0].imag() * a[1].imag()) / mo;
cx[0] = complex<_Ty>(sb, xb);
return(1);
}
pq = Abs(a[k]);
}
q = log(pq);
q = q / (1.0 * k);
q = exp(q);
p = q;
w = w * p;
for(i = 1; i <= k; i++)
{
a[i] /= q;
q *= p;
}
x = 0.0001;
x1 = x;
y = 0.2;
y1 = y;
dx = 1.0;
g = 1.0e+37;
zjn:
u = a[0].real();
v = a[0].imag();
for(i = 1; i <= k; i ++)
{
p = u * x1;
q = v * y1;
pq = (u + v) * (x1 + y1);
u = p - q + a[i].real();
v = pq - p - q + a[i].imag();
}
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 zjn;
}
else
{
gg60(t, x, y, x1, y1, dx, dy, p, q, k, it);
if(t >= 1.0e-03) goto zjn;
if(g > 1.0e-18)
{
it = 0;
g65(x, y, x1, y1, dx, dy, dd, dc, c, k, is, it);
if(it == 0) goto zjn;
}
}
gg90(cx, a, x, y, p, q, w, k);
}
else
{
g = g1;
x = x1;
y = y1;
is = 0;
if(g <= 1.0e-22)
gg90(cx, a, x, y, p, q, w, k);
else
{
u1 = k * a[0].real();
v1 = a[0].imag();
for(i = 2; i <= k; i++)
{
p = u1 * x;
q = v1 * y;
pq = (u1 + v1) * (x + y);
u1 = p - q + (k - i + 1) * a[i - 1].real();
v1 = pq - p - q + (k - i + 1) * a[i - 1].imag();
}
p = u1 * u1 + v1 * v1;
if(p <= 1.0e-20)
{
it = 0;
g65(x, y, x1, y1, dx, dy, dd, dc, c, k, is, it);
if(it == 0) goto zjn;
gg90(cx, a, x, y, p, q, w, k);
}
else
{
dx = (u * u1 + v * v1) / p;
dy = (u1 * v - v1 * u) / p;
t = 1.0 + 4.0 / k;
gg60(t, x, y, x1, y1, dx, dy, p, q, k, it);
if(t >= 1.0e-03) goto zjn;
if(g > 1.0e-18)
?? 快捷鍵說明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -