开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 数学计算 > 查看源码
NLequation.cpp
资源名称：科学与工程数值计算算法配套源码vc++.zip　[点击查看]
上传用户：weigute
上传日期：2007-03-02
资源大小：1287k
文件大小：31k
源码类别：

数学计算

开发平台：

Visual C++

NLequation.cpp：源码内容
  1. //////////////////////////////////////////////////////////////////////
  2. // NLequation.cpp
  3. //
  4. // 求解线性方程组的类 CNLequation 的实现代码
  5. //
  6. // 周长发编制, 2002/8
  7. //////////////////////////////////////////////////////////////////////
  9. #include "stdafx.h"
  10. #include "NLequation.h"
  11. #include "LEquations.h"
  13. #ifdef _DEBUG
  14. #undef THIS_FILE
  15. static char THIS_FILE[]=__FILE__;
  16. #define new DEBUG_NEW
  17. #endif
  19. //////////////////////////////////////////////////////////////////////
  20. // Construction/Destruction
  21. //////////////////////////////////////////////////////////////////////
  23. //////////////////////////////////////////////////////////////////////
  24. // 基本构造函数
  25. //////////////////////////////////////////////////////////////////////
  26. CNLequation::CNLequation()
  27. {
  28. }
  30. //////////////////////////////////////////////////////////////////////
  31. // 析构函数
  32. //////////////////////////////////////////////////////////////////////
  33. CNLequation::~CNLequation()
  34. {
  35. }
  37. //////////////////////////////////////////////////////////////////////
  38. // 求非线性方程实根的对分法
  39. //
  40. // 调用时，须覆盖计算方程左端函数f(x)值的虚函数double Func(double x)
  41. //
  42. // 参数：
  43. // 1. int nNumRoots - 在[xStart, xEnd]内实根个数的预估值
  44. // 2. double x[] - 一维数组，长度为m。返回在区间[xStart, xEnd]内搜索到
  45. //                 的实根，实根个数由函数值返回
  46. // 3. double xStart - 求根区间的左端点
  47. // 4. double xEnd - 求根区间的右端点
  48. // 5. double dblStep - 搜索求根时采用的步长
  49. // 6. double eps - 精度控制参数，默认值为0.000001
  50. //
  51. // 返回值：int 型，求得的实根的数目
  52. //////////////////////////////////////////////////////////////////////
  53. int CNLequation::GetRootBisect(int nNumRoots, double x[], double xStart, double xEnd, double dblStep, double eps /*= 0.000001*/)
  54. {
  55. int n,js;
  56.     double z,y,z1,y1,z0,y0;
  58. // 根的个数清0
  59.     n = 0; 
  61. // 从左端点开始搜索
  62. z = xStart; 
  63. y = Func(z);
  65. // 循环求解
  66.     while ((z<=xEnd+dblStep/2.0)&&(n!=nNumRoots))
  67.     { 
  68. if (fabs(y)<eps)
  69.         { 
  70. n=n+1; 
  71. x[n-1]=z;
  72.             z=z+dblStep/2.0; 
  73. y=Func(z);
  74.         }
  75.         else
  76.         { 
  77. z1=z+dblStep; 
  78. y1=Func(z1);
  79.             
  80. if (fabs(y1)<eps)
  81.             { 
  82. n=n+1; 
  83. x[n-1]=z1;
  84.                 z=z1+dblStep/2.0; 
  85. y=Func(z);
  86.             }
  87.             else if (y*y1>0.0)
  88.             { 
  89. y=y1; 
  90. z=z1;
  91. }
  92.             else
  93.             { 
  94. js=0;
  95.                 while (js==0)
  96.                 { 
  97. if (fabs(z1-z)<eps)
  98.                     { 
  99. n=n+1; 
  100. x[n-1]=(z1+z)/2.0;
  101.                         z=z1+dblStep/2.0; y=Func(z);
  102.                         js=1;
  103.                     }
  104.                     else
  105.                     { 
  106. z0=(z1+z)/2.0; 
  107. y0=Func(z0);
  108.                         if (fabs(y0)<eps)
  109.                         { 
  110. x[n]=z0; 
  111. n=n+1; 
  112. js=1;
  113.                             z=z0+dblStep/2.0; 
  114. y=Func(z);
  115.                         }
  116.                         else if ((y*y0)<0.0)
  117.                         { 
  118. z1=z0; 
  119. y1=y0;
  120. }
  121.                         else 
  122. { 
  123. z=z0; 
  124. y=y0;
  125. }
  126.                     }
  127.                 }
  128.             }
  129.         }
  130.     }
  131.     
  132. // 返回实根的数目
  133. return(n);
  134. }
  136. //////////////////////////////////////////////////////////////////////
  137. // 求非线性方程一个实根的牛顿法
  138. //
  139. // 调用时，须覆盖计算方程左端函数f(x)及其一阶导数f'(x)值的虚函数
  140. //                  void Func(double x, double y[])
  141. //         y(0) 返回f(x)的值
  142. //         y(1) 返回f'(x)的值
  143. //
  144. // 参数：
  145. // 1. double *x - 传入迭代初值（猜测解），返回在区间求得的一个实根
  146. // 2. int nMaxIt - 递归次数，默认值为60
  147. // 3. double eps - 精度控制参数，默认值为0.000001
  148. //
  149. // 返回值：BOOL 型，求解是否成功
  150. //////////////////////////////////////////////////////////////////////
  151. BOOL CNLequation::GetRootNewton(double* x, int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
  152. { 
  153.     int l;
  154.     double y[2],d,p,x0,x1;
  156. // 条件值
  157.     l=nMaxIt; 
  158. x0=*x;
  159.     Func(x0,y);
  160.     
  161. // 求解，控制精度
  162. d=eps+1.0;
  163.     while ((d>=eps)&&(l!=0))
  164.     { 
  165. if (y[1] == 0.0)
  166. return FALSE;
  168.         x1=x0-y[0]/y[1];
  169.         Func(x1,y);
  170.         
  171. d=fabs(x1-x0); 
  172. p=fabs(y[0]);
  173.         if (p>d) 
  174. d=p;
  175.         x0=x1; 
  176. l=l-1;
  177.     }
  178.     
  179. *x=x1;
  181. return TRUE;
  182. }
  184. //////////////////////////////////////////////////////////////////////
  185. // 求非线性方程一个实根的埃特金迭代法
  186. //
  187. // 调用时，须覆盖计算方程左端函数f(x)值的虚函数double Func(double x)
  188. //
  189. // 参数：
  190. // 1. double *x - 传入迭代初值（猜测解），返回在区间求得的一个实根
  191. // 2. int nMaxIt - 递归次数，默认值为60
  192. // 3. double eps - 精度控制参数，默认值为0.000001
  193. //
  194. // 返回值：BOOL 型，求解是否成功
  195. //////////////////////////////////////////////////////////////////////
  196. BOOL CNLequation::GetRootAitken(double* x, int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
  197. { 
  198. int flag,l;
  199.     double u,v,x0;
  201. // 求解条件
  202.     l=0; 
  203. x0=*x; 
  204. flag=0;
  206. // 迭代求解
  207.     while ((flag==0)&&(l!=nMaxIt))
  208.     { 
  209. l=l+1; 
  210.         u=Func(x0); 
  211. v=Func(u);
  212.         if (fabs(u-v)<eps) 
  213. { 
  214. x0=v; 
  215. flag=1; 
  216. }
  217.         else 
  218. x0=v-(v-u)*(v-u)/(v-2.0*u+x0);
  219.     }
  220.     
  221. *x=x0; 
  222.     
  223. // 是否在指定的迭代次数内达到求解精度
  224. return (nMaxIt > l);
  225. }
  227. //////////////////////////////////////////////////////////////////////
  228. // 求非线性方程一个实根的连分式解法
  229. //
  230. // 调用时，须覆盖计算方程左端函数f(x)值的虚函数double Func(double x)
  231. //
  232. // 参数：
  233. // 1. double *x - 传入迭代初值（猜测解），返回在区间求得的一个实根
  234. // 2. double eps - 精度控制参数，默认值为0.000001
  235. //
  236. // 返回值：BOOL 型，求解是否成功
  237. //////////////////////////////////////////////////////////////////////
  238. BOOL CNLequation::GetRootPq(double* x, double eps /*= 0.000001*/)
  239. { 
  240.     int i,j,m,it,l;
  241.     double a[10],y[10],z,h,x0,q;
  243. // 求解条件
  244.     l=10; 
  245. q=1.0e+35; 
  246. x0=*x; 
  247. h=0.0;
  248.     
  249. // 连分式求解
  250. while (l!=0)
  251.     { 
  252. l=l-1; 
  253. j=0; 
  254. it=l;
  255.         while (j<=7)
  256.         { 
  257. if (j<=2) 
  258. z=x0+0.1*j;
  259.             else 
  260. z=h;
  262. y[j]=Func(z);
  263.             h=z;
  264.             if (j==0) 
  265. a[0]=z;
  266.             else
  267.             { 
  268. m=0; 
  269. i=0;
  270.                 while ((m==0)&&(i<=j-1))
  271.                 { 
  272. if (fabs(h-a[i])+1.0==1.0) 
  273. m=1;
  274.                     else 
  275. h=(y[j]-y[i])/(h-a[i]);
  276.                      
  277. i=i+1;
  278.                 }
  279.                 a[j]=h;
  280.                 if (m!=0) 
  281. a[j]=q;
  282.                 h=0.0;
  283.                 for (i=j-1; i>=0; i--)
  284.                 { 
  285. if (fabs(a[i+1]+h)+1.0==1.0) 
  286. h=q;
  287.                     else 
  288. h=-y[i]/(a[i+1]+h);
  289.                 }
  290.                  
  291. h=h+a[0];
  292.             }
  293.              
  294. if (fabs(y[j])>=eps) 
  295. j=j+1;
  296.             else 
  297. { 
  298. j=10; 
  299. l=0;
  300. }
  301. }
  302.         
  303. x0=h;
  304. }
  306.     *x=h;
  307.     
  308. // 是否在10阶连分式内求的实根？
  309. return (10>it);
  310. }
  312. //////////////////////////////////////////////////////////////////////
  313. // 求实系数代数方程全部根的QR方法
  314. //
  315. // 参数：
  316. // 1. int n - 多项式方程的次数
  317. // 2. double dblCoef[] - 一维数组，长度为n+1，按降幂次序依次存放n次多项式方程的n+1个系数
  318. // 3. double xr[] - 一维数组，长度为n，返回n个根的实部
  319. // 4. double xi[] - 一维数组，长度为n，返回n个根的虚部
  320. // 5. int nMaxIt - 迭代次数，默认值为 60
  321. // 6. double eps - 精度控制参数，默认值为 0.000001
  322. //
  323. // 返回值：BOOL 型，求解是否成功
  324. //////////////////////////////////////////////////////////////////////
  325. BOOL CNLequation::GetRootQr(int n, double dblCoef[], double xr[], double xi[], int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
  326. { 
  327. // 初始化矩阵
  328. CMatrix mtxQ;
  329. mtxQ.Init(n, n);
  330.     double *q = mtxQ.GetData();
  332. // 构造赫申伯格矩阵
  333.     for (int j=0; j<=n-1; j++)
  334. q[j]=-dblCoef[n-j-1]/dblCoef[n];
  336.     for (j=n; j<=n*n-1; j++)
  337. q[j]=0.0;
  339.     for (int i=0; i<=n-2; i++)
  340. q[(i+1)*n+i]=1.0;
  342. // 求赫申伯格矩阵的特征值和特征向量，即为方程的解
  343.     if (mtxQ.HBergEigenv(xr, xi, nMaxIt, eps))
  344. return TRUE;
  346. return FALSE;
  347. }
  349. //////////////////////////////////////////////////////////////////////
  350. // 求实系数代数方程全部根的牛顿下山法
  351. //
  352. // 参数：
  353. // 1. int n - 多项式方程的次数
  354. // 2. double dblCoef[] - 一维数组，长度为n+1，按降幂次序依次存放n次多项式方程的n+1个系数
  355. // 3. double xr[] - 一维数组，长度为n，返回n个根的实部
  356. // 4. double xi[] - 一维数组，长度为n，返回n个根的虚部
  357. //
  358. // 返回值：BOOL 型，求解是否成功
  359. //////////////////////////////////////////////////////////////////////
  360. BOOL CNLequation::GetRootNewtonDownHill(int n, double dblCoef[], double xr[], double xi[])
  361. { 
  362. int m,i,jt,k,is,it;
  363.     double t,x,y,x1,y1,dx,dy,p,q,w,dd,dc,c;
  364.     double g,u,v,pq,g1,u1,v1;
  365.     
  366. // 初始判断
  367.     m=n;
  368.     while ((m>0)&&(fabs(dblCoef[m])+1.0==1.0)) 
  369. m=m-1;
  371. // 求解失败
  372.     if (m<=0)
  373. return FALSE;
  375.     for (i=0; i<=m; i++)
  376. dblCoef[i]=dblCoef[i]/dblCoef[m];
  377.     
  378. for (i=0; i<=m/2; i++)
  379.     { 
  380. w=dblCoef[i]; 
  381. dblCoef[i]=dblCoef[m-i]; 
  382. dblCoef[m-i]=w;
  383. }
  384.     
  385. // 迭代求解
  386. k=m; 
  387. is=0; 
  388. w=1.0;
  389.     jt=1;
  390.     while (jt==1)
  391.     { 
  392. pq=fabs(dblCoef[k]);
  393. while (pq<1.0e-12)
  394.         { 
  395. xr[k-1]=0.0; 
  396. xi[k-1]=0.0; 
  397. k=k-1;
  398.             if (k==1)
  399.             { 
  400. xr[0]=-dblCoef[1]*w/dblCoef[0]; 
  401. xi[0]=0.0;
  402.                 
  403. return TRUE;
  404.             }
  405.             
  406. pq=fabs(dblCoef[k]);
  407.         }
  409. q=log(pq); 
  410. q=q/(1.0*k); 
  411. q=exp(q);
  412.         p=q; 
  413. w=w*p;
  414.         for (i=1; i<=k; i++)
  415.         { 
  416. dblCoef[i]=dblCoef[i]/q; 
  417. q=q*p;
  418. }
  419.         
  420. x=0.0001; 
  421. x1=x; 
  422. y=0.2; 
  423. y1=y; 
  424. dx=1.0;
  425.         g=1.0e+37; 
  426. l40:
  427.         u=dblCoef[0]; v=0.0;
  428.         for (i=1; i<=k; i++)
  429.         { 
  430. p=u*x1; 
  431. q=v*y1;
  432.             pq=(u+v)*(x1+y1);
  433.             u=p-q+dblCoef[i]; 
  434. v=pq-p-q;
  435.         }
  436.         
  437. g1=u*u+v*v;
  438.         if (g1>=g)
  439.         { 
  440. if (is!=0)
  441.             { 
  442. it=1;
  443.                 g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
  444.                 if (it==0) 
  445. goto l40;
  446.             }
  447.             else
  448.             { 
  449. g60(&t,&x,&y,&x1,&y1,&dx,&dy,&p,&q,&k,&it);
  450.                 if (t>=1.0e-03) 
  451. goto l40;
  452.                 
  453. if (g>1.0e-18)
  454.                 { 
  455. it=0;
  456.                     g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
  457.                     if (it==0) 
  458. goto l40;
  459.                 }
  460.             }
  461.             
  462. g90(xr,xi,dblCoef,&x,&y,&p,&q,&w,&k);
  463.         }
  464.         else
  465.         { 
  466. g=g1; 
  467. x=x1; 
  468. y=y1; 
  469. is=0;
  470.             if (g<=1.0e-22)
  471. g90(xr,xi,dblCoef,&x,&y,&p,&q,&w,&k);
  472.             else
  473.             { 
  474. u1=k*dblCoef[0]; 
  475. v1=0.0;
  476.                 for (i=2; i<=k; i++)
  477.                 { 
  478. p=u1*x; 
  479. q=v1*y; 
  480. pq=(u1+v1)*(x+y);
  481.                     u1=p-q+(k-i+1)*dblCoef[i-1];
  482.                     v1=pq-p-q;
  483.                 }
  484.                 
  485. p=u1*u1+v1*v1;
  486.                 if (p<=1.0e-20)
  487.                 { 
  488. it=0;
  489.                     g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
  490.                     if (it==0) 
  491. goto l40;
  493.                     g90(xr,xi,dblCoef,&x,&y,&p,&q,&w,&k);
  494.                 }
  495.                 else
  496.                 { 
  497. dx=(u*u1+v*v1)/p;
  498.                     dy=(u1*v-v1*u)/p;
  499.                     t=1.0+4.0/k;
  500.                     g60(&t,&x,&y,&x1,&y1,&dx,&dy,&p,&q,&k,&it);
  501.                     if (t>=1.0e-03) 
  502. goto l40;
  504.                     if (g>1.0e-18)
  505.                     { 
  506. it=0;
  507.                         g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
  508.                         if (it==0) 
  509. goto l40;
  510.                     }
  511.                     
  512. g90(xr,xi,dblCoef,&x,&y,&p,&q,&w,&k);
  513.                 }
  514.             }
  515.         }
  516.         
  517. if (k==1) 
  518. jt=0;
  519.         else 
  520. jt=1;
  521.     }
  522.     
  523. return TRUE;
  524. }
  527. //////////////////////////////////////////////////////////////////////
  528. // 内部函数
  529. //////////////////////////////////////////////////////////////////////
  530. void CNLequation::g60(double* t,double* x,double* y,double* x1,double* y1,double* dx,double* dy,double* p,double* q,int* k,int* it)
  531. { 
  532. *it=1;
  533.     while (*it==1)
  534.     { 
  535. *t=*t/1.67; 
  536. *it=0;
  537.         *x1=*x-(*t)*(*dx);
  538.         *y1=*y-(*t)*(*dy);
  539.         if (*k>=50)
  540. { 
  541. *p=sqrt((*x1)*(*x1)+(*y1)*(*y1));
  542.             *q=exp(85.0/(*k));
  543.             if (*p>=*q) 
  544. *it=1;
  545.         }
  546.     }
  547. }
  549. //////////////////////////////////////////////////////////////////////
  550. // 内部函数
  551. //////////////////////////////////////////////////////////////////////
  552. void CNLequation::g90(double xr[],double xi[],double dblCoef[],double* x,double* y,double* p,double* q,double* w,int* k)
  553. { 
  554. int i;
  555.     
  556. if (fabs(*y)<=1.0e-06)
  557.     { 
  558. *p=-(*x); 
  559. *y=0.0; 
  560. *q=0.0;
  561. }
  562.     else
  563.     { 
  564. *p=-2.0*(*x); 
  565. *q=(*x)*(*x)+(*y)*(*y);
  566.         xr[*k-1]=(*x)*(*w);
  567.         xi[*k-1]=-(*y)*(*w);
  568.         *k=*k-1;
  569.     }
  570.     
  571. for (i=1; i<=*k; i++)
  572.     { 
  573. dblCoef[i]=dblCoef[i]-dblCoef[i-1]*(*p);
  574.         dblCoef[i+1]=dblCoef[i+1]-dblCoef[i-1]*(*q);
  575.     }
  576.     
  577. xr[*k-1]=(*x)*(*w); 
  578. xi[*k-1]=(*y)*(*w);
  579.     *k=*k-1;
  580.     if (*k==1)
  581.     { 
  582. xr[0]=-dblCoef[1]*(*w)/dblCoef[0]; 
  583. xi[0]=0.0;
  584. }
  585. }
  587. //////////////////////////////////////////////////////////////////////
  588. // 内部函数
  589. //////////////////////////////////////////////////////////////////////
  590. void CNLequation::g65(double* x,double* y,double* x1,double* y1,double* dx,double* dy,double* dd,double* dc,double* c,int* k,int* is,int* it)
  591. { 
  592. if (*it==0)
  593.     { 
  594. *is=1;
  595.         *dd=sqrt((*dx)*(*dx)+(*dy)*(*dy));
  596.         if (*dd>1.0) 
  597. *dd=1.0;
  598.         *dc=6.28/(4.5*(*k)); 
  599. *c=0.0;
  600.     }
  601.     
  602. while(TRUE)
  603.     { 
  604. *c=*c+(*dc);
  605.         *dx=(*dd)*cos(*c); 
  606. *dy=(*dd)*sin(*c);
  607.         *x1=*x+*dx; 
  608. *y1=*y+*dy;
  609.         if (*c<=6.29)
  610.         { 
  611. *it=0; 
  612. return;
  613. }
  614.         
  615. *dd=*dd/1.67;
  616.         if (*dd<=1.0e-07)
  617.         { 
  618. *it=1; 
  619. return;
  620. }
  621.         
  622. *c=0.0;
  623.     }
  624. }
  626. //////////////////////////////////////////////////////////////////////
  627. // 求复系数代数方程全部根的牛顿下山法
  628. //
  629. // 参数：
  630. // 1. int n - 多项式方程的次数
  631. // 2. double ar[] - 一维数组，长度为n+1，按降幂次序依次存放n次多项式方程的n+1个系数的实部
  632. // 3. double ai[] - 一维数组，长度为n+1，按降幂次序依次存放n次多项式方程的n+1个系数的虚部
  633. // 4. double xr[] - 一维数组，长度为n，返回n个根的实部
  634. // 5. double xi[] - 一维数组，长度为n，返回n个根的虚部
  635. //
  636. // 返回值：BOOL 型，求解是否成功
  637. //////////////////////////////////////////////////////////////////////
  638. BOOL CNLequation::GetRootNewtonDownHill(int n, double ar[], double ai[], double xr[], double xi[])
  639. { 
  640. int m,i,jt,k,is,it;
  641.     double t,x,y,x1,y1,dx,dy,p,q,w,dd,dc,c;
  642.     double g,u,v,pq,g1,u1,v1;
  644. // 初始判断
  645.     m=n;
  646.     p=sqrt(ar[m]*ar[m]+ai[m]*ai[m]);
  647.     while ((m>0)&&(p+1.0==1.0))
  648.     {  
  649. m=m-1;
  650.         p=sqrt(ar[m]*ar[m]+ai[m]*ai[m]);
  651.     }
  652.     
  653. // 求解失败
  654. if (m<=0)
  655. return FALSE;
  657.     for (i=0; i<=m; i++)
  658.     { 
  659. ar[i]=ar[i]/p; 
  660. ai[i]=ai[i]/p;
  661. }
  662.     
  663. for (i=0; i<=m/2; i++)
  664.     { 
  665. w=ar[i]; 
  666. ar[i]=ar[m-i]; 
  667. ar[m-i]=w;
  668.         w=ai[i]; 
  669. ai[i]=ai[m-i]; 
  670. ai[m-i]=w;
  671.     }
  672.     
  673. // 迭代求解
  674. k=m; 
  675. is=0; 
  676. w=1.0;
  677.     jt=1;
  678.     while (jt==1)
  679.     { 
  680. pq=sqrt(ar[k]*ar[k]+ai[k]*ai[k]);
  681. while (pq<1.0e-12)
  682.         { 
  683. xr[k-1]=0.0; 
  684. xi[k-1]=0.0; 
  685. k=k-1;
  686.             if (k==1)
  687.             { 
  688. p=ar[0]*ar[0]+ai[0]*ai[0];
  689.                 xr[0]=-w*(ar[0]*ar[1]+ai[0]*ai[1])/p;
  690.                 xi[0]=w*(ar[1]*ai[0]-ar[0]*ai[1])/p;
  691.                 
  692. return TRUE;
  693.             }
  694.             
  695. pq=sqrt(ar[k]*ar[k]+ai[k]*ai[k]);
  696.         }
  698. q=log(pq); 
  699. q=q/(1.0*k); 
  700. q=exp(q);
  701.         p=q; 
  702. w=w*p;
  703.         for (i=1; i<=k; i++)
  704.         { 
  705. ar[i]=ar[i]/q; 
  706. ai[i]=ai[i]/q; 
  707. q=q*p;
  708. }
  709.         
  710. x=0.0001; 
  711. x1=x; 
  712. y=0.2; 
  713. y1=y; 
  714. dx=1.0;
  715.         g=1.0e+37; 
  716. l40:
  717.         u=ar[0]; 
  718. v=ai[0];
  719.         for (i=1; i<=k; i++)
  720.         { 
  721. p=u*x1; 
  722. q=v*y1;
  723.             pq=(u+v)*(x1+y1);
  724.             u=p-q+ar[i]; 
  725. v=pq-p-q+ai[i];
  726.         }
  727.         
  728. g1=u*u+v*v;
  729.         if (g1>=g)
  730.         { 
  731. if (is!=0)
  732.             { 
  733. it=1;
  734.                 g65c(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
  735.                 if (it==0) 
  736. goto l40;
  737.             }
  738.             else
  739.             { 
  740. g60c(&t,&x,&y,&x1,&y1,&dx,&dy,&p,&q,&k,&it);
  741.                 if (t>=1.0e-03) 
  742. goto l40;
  743.                 
  744. if (g>1.0e-18)
  745.                 { 
  746. it=0;
  747.                     g65c(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
  748.                     if (it==0) 
  749. goto l40;
  750.                 }
  751.             }
  752.             
  753. g90c(xr,xi,ar,ai,&x,&y,&p,&w,&k);
  754.         }
  755.         else
  756.         { 
  757. g=g1; 
  758. x=x1; 
  759. y=y1; 
  760. is=0;
  761.             if (g<=1.0e-22)
  762. g90c(xr,xi,ar,ai,&x,&y,&p,&w,&k);
  763.             else
  764.             { 
  765. u1=k*ar[0]; 
  766. v1=ai[0];
  767.                 for (i=2; i<=k; i++)
  768.                 { 
  769. p=u1*x; 
  770. q=v1*y; 
  771. pq=(u1+v1)*(x+y);
  772.                     u1=p-q+(k-i+1)*ar[i-1];
  773.                     v1=pq-p-q+(k-i+1)*ai[i-1];
  774.                 }
  775.                 
  776. p=u1*u1+v1*v1;
  777.                 if (p<=1.0e-20)
  778.                 { 
  779. it=0;
  780.                     g65c(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
  781.                     if (it==0) 
  782. goto l40;
  783.                     
  784. g90c(xr,xi,ar,ai,&x,&y,&p,&w,&k);
  785.                 }
  786.                 else
  787.                 { 
  788. dx=(u*u1+v*v1)/p;
  789.                     dy=(u1*v-v1*u)/p;
  790.                     t=1.0+4.0/k;
  791.                     g60c(&t,&x,&y,&x1,&y1,&dx,&dy,&p,&q,&k,&it);
  792.                     if (t>=1.0e-03) 
  793. goto l40;
  794.                     
  795. if (g>1.0e-18)
  796.                     { 
  797. it=0;
  798.                         g65c(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
  799.                         if (it==0) 
  800. goto l40;
  801.                     }
  802.                     
  803. g90c(xr,xi,ar,ai,&x,&y,&p,&w,&k);
  804.                 }
  805.             }
  806.         }
  807.         
  808. if (k==1) 
  809. jt=0;
  810.         else 
  811. jt=1;
  812.     }
  813.     
  814. return TRUE;
  815. }
  817. //////////////////////////////////////////////////////////////////////
  818. // 内部函数
  819. //////////////////////////////////////////////////////////////////////
  820. void CNLequation::g60c(double* t,double* x,double* y,double* x1,double* y1,double* dx,double* dy,double* p,double* q,int* k,int* it)
  821. { 
  822. *it=1;
  823.     while (*it==1)
  824.     { 
  825. *t=*t/1.67; 
  826. *it=0;
  827.         *x1=*x-*t*(*dx);
  828.         *y1=*y-*t*(*dy);
  829.         if (*k>=30)
  830. { 
  831. *p=sqrt(*x1*(*x1)+*y1*(*y1));
  832.             *q=exp(75.0/(*k));
  833.             if (*p>=*q) 
  834. *it=1;
  835.         }
  836.     }
  837. }
  839. //////////////////////////////////////////////////////////////////////
  840. // 内部函数
  841. //////////////////////////////////////////////////////////////////////
  842. void CNLequation::g90c(double xr[],double xi[],double ar[],double ai[],double* x,double* y,double* p,double* w,int* k)
  843. { 
  844. int i;
  845.     for (i=1; i<=*k; i++)
  846.     { 
  847. ar[i]=ar[i]+ar[i-1]*(*x)-ai[i-1]*(*y);
  848.         ai[i]=ai[i]+ar[i-1]*(*y)+ai[i-1]*(*x);
  849.     }
  850.     
  851. xr[*k-1]=*x*(*w); 
  852. xi[*k-1]=*y*(*w);
  853.     *k=*k-1;
  854.     if (*k==1)
  855.     { 
  856. *p=ar[0]*ar[0]+ai[0]*ai[0];
  857.         xr[0]=-*w*(ar[0]*ar[1]+ai[0]*ai[1])/(*p);
  858.         xi[0]=*w*(ar[1]*ai[0]-ar[0]*ai[1])/(*p);
  859.     }
  860. }
  862. //////////////////////////////////////////////////////////////////////
  863. // 内部函数
  864. //////////////////////////////////////////////////////////////////////
  865. void CNLequation::g65c(double* x,double* y,double* x1,double* y1,double* dx,double* dy,double* dd,double* dc,double* c,int* k,int* is,int* it)
  866. { 
  867. if (*it==0)
  868.     { 
  869. *is=1;
  870.         *dd=sqrt(*dx*(*dx)+*dy*(*dy));
  871.         if (*dd>1.0) 
  872. *dd=1.0;
  873.         *dc=6.28/(4.5*(*k)); 
  874. *c=0.0;
  875.     }
  876.     
  877. while(TRUE)
  878.     { 
  879. *c=*c+*dc;
  880.         *dx=*dd*cos(*c); 
  881. *dy=*dd*sin(*c);
  882.         *x1=*x+*dx; 
  883. *y1=*y+*dy;
  884.         if (*c<=6.29)
  885.         { 
  886. *it=0; 
  887. return;
  888. }
  889.         
  890. *dd=*dd/1.67;
  891.         if (*dd<=1.0e-07)
  892.         { 
  893. *it=1; 
  894. return;
  895. }
  896.         
  897. *c=0.0;
  898.     }
  899. }
  901. //////////////////////////////////////////////////////////////////////
  902. // 求非线性方程一个实根的蒙特卡洛法
  903. //
  904. // 调用时，须覆盖计算方程左端函数f(x)值的虚函数double Func(double x)
  905. //
  906. // 参数：
  907. // 1. double* x - 传入初值（猜测解），返回求得的实根
  908. // 2. double xStart - 均匀分布的端点初值
  909. // 3. int nControlB - 控制参数
  910. // 4. double eps - 控制精度，默认值为0.000001
  911. //
  912. // 返回值：无
  913. //////////////////////////////////////////////////////////////////////
  914. void CNLequation::GetRootMonteCarlo(double* x, double xStart, int nControlB, double eps /*= 0.000001*/)
  915. { 
  916. int k;
  917.     double xx,a,y,x1,y1,r;
  918.     
  919. // 求解条件
  920. a = xStart; 
  921. k = 1; 
  922. r = 1.0;
  924. // 初值
  925. xx = *x; 
  926. y = Func(xx);
  928. // 精度控制求解
  929.     while (a>=eps)
  930.     { 
  931. x1=rnd(&r);
  933. x1=-a+2.0*a*x1;
  934.         x1=xx+x1; 
  935. y1=Func(x1);
  936.         
  937. k=k+1;
  938.         if (fabs(y1)>=fabs(y))
  939.         { 
  940. if (k>nControlB) 
  941. { 
  942. k=1; 
  943. a=a/2.0; 
  944. }
  945. }
  946.         else
  947.         { 
  948. k=1; 
  949. xx=x1; 
  950. y=y1;
  951.             if (fabs(y)<eps)
  952.             { 
  953. *x=xx; 
  954. return; 
  955. }
  956.         }
  957.     }
  958.     
  959. *x=xx;
  960. }
  962. //////////////////////////////////////////////////////////////////////
  963. // 内部函数
  964. //////////////////////////////////////////////////////////////////////
  965. double CNLequation::rnd(double* r)
  966. {
  967. int m;
  968.     double s,u,v,p;
  969.     
  970. s=65536.0; 
  971. u=2053.0; 
  972. v=13849.0;
  973.     m=(int)(*r/s); 
  974. *r=*r-m*s;
  975.     *r=u*(*r)+v; 
  976. m=(int)(*r/s);
  977.     *r=*r-m*s; p=*r/s;
  978.     
  979. return(p);
  980. }
  982. //////////////////////////////////////////////////////////////////////
  983. // 求实函数或复函数方程一个复根的蒙特卡洛法
  984. //
  985. // 调用时，须覆盖计算方程左端函数的模值||f(x, y)||的虚函数
  986. //           double Func(double x, double y)
  987. //
  988. // 参数：
  989. // 1. double* x - 传入初值（猜测解）的实部，返回求得的根的实部
  990. // 2. double* y - 传入初值（猜测解）的虚部，返回求得的根的虚部
  991. // 3. double xStart - 均匀分布的端点初值
  992. // 4. int nControlB - 控制参数
  993. // 5. double eps - 控制精度，默认值为0.000001
  994. // 6. double xi[] - 一维数组，长度为n，返回n个根的虚部
  995. //
  996. // 返回值：无
  997. //////////////////////////////////////////////////////////////////////
  998. void CNLequation::GetRootMonteCarlo(double* x, double* y, double xStart, int nControlB, double eps /*= 0.000001*/)
  999. { 
  1000.     int k;
  1001.     double xx,yy,a,r,z,x1,y1,z1;
  1003. // 求解条件与初值
  1004.     a=xStart; 
  1005. k=1; 
  1006. r=1.0; 
  1007. xx=*x; 
  1008. yy=*y;
  1009.     z=Func(xx,yy);
  1010.     
  1011. // 精度控制求解
  1012. while (a>=eps)
  1013.     { 
  1014. x1=-a+2.0*a*rnd(&r); 
  1015. x1=xx+x1; 
  1016.         y1=-a+2.0*a*rnd(&r); 
  1017. y1=yy+y1;
  1018.         
  1019. z1=Func(x1,y1);
  1020.         
  1021. k=k+1;
  1022.         if (z1>=z)
  1023.         { 
  1024. if (k>nControlB) 
  1025. { 
  1026. k=1; 
  1027. a=a/2.0; 
  1028. }
  1029. }
  1030.         else
  1031.         { 
  1032. k=1; 
  1033. xx=x1; 
  1034. yy=y1;  
  1035. z=z1;
  1036.             if (z<eps)
  1037.             { 
  1038. *x=xx; 
  1039. *y=yy; 
  1040. return; 
  1041. }
  1042.         }
  1043.     }
  1044.     
  1045. *x=xx; 
  1046. *y=yy; 
  1047. }
  1049. //////////////////////////////////////////////////////////////////////
  1050. // 求非线性方程组一组实根的梯度法
  1051. //
  1052. // 调用时，须覆盖计算方程左端函数f(x)值及其偏导数值的虚函数
  1053. //         double Func(double x[], double y[])
  1054. //
  1055. // 参数：
  1056. // 1. int n - 方程的个数，也是未知数的个数
  1057. // 2. double x[] - 一维数组，长度为n，存放一组初值x0, x1, …, xn-1，
  1058. //                 返回时存放方程组的一组实根
  1059. // 3. int nMaxIt - 迭代次数，默认值为60
  1060. // 4. double eps - 控制精度，默认值为0.000001
  1061. //
  1062. // 返回值：BOOL 型，求解是否成功
  1063. //////////////////////////////////////////////////////////////////////
  1064. BOOL CNLequation::GetRootsetGrad(int n, double x[], int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
  1065. { 
  1066.     int l,j;
  1067.     double f,d,s,*y;
  1069.     y = new double[n];
  1071.     l=nMaxIt;
  1072.     f=Func(x,y);
  1074. // 控制精度，迭代求解
  1075.     while (f>=eps)
  1076.     { 
  1077. l=l-1;
  1078.         if (l==0) 
  1079. { 
  1080. delete[] y; 
  1081. return TRUE;
  1082. }
  1083.         
  1084. d=0.0;
  1085.         for (j=0; j<=n-1; j++) 
  1086. d=d+y[j]*y[j];
  1087.         if (d+1.0==1.0) 
  1088. { 
  1089. delete[] y; 
  1090. return FALSE;
  1091. }
  1092.         
  1093. s=f/d;
  1094.         for (j=0; j<=n-1; j++) 
  1095. x[j]=x[j]-s*y[j];
  1096.         
  1097. f=Func(x,y);
  1098.     }
  1099.     
  1100. delete[] y; 
  1102. // 是否在有效迭代次数内达到精度
  1103. return (nMaxIt>l);
  1104. }
  1106. //////////////////////////////////////////////////////////////////////
  1107. // 求非线性方程组一组实根的拟牛顿法
  1108. //
  1109. // 调用时，须覆盖计算方程左端函数f(x)值及其偏导数值的虚函数
  1110. //         double Func(double x[], double y[])
  1111. //
  1112. // 参数：
  1113. // 1. int n - 方程的个数，也是未知数的个数
  1114. // 2. double x[] - 一维数组，长度为n，存放一组初值x0, x1, …, xn-1，
  1115. //                 返回时存放方程组的一组实根
  1116. // 3. double t - 控制h大小的变量，0<t<1
  1117. // 4. double h - 增量初值
  1118. // 3. int nMaxIt - 迭代次数，默认值为60
  1119. // 4. double eps - 控制精度，默认值为0.000001
  1120. //
  1121. // 返回值：BOOL 型，求解是否成功
  1122. //////////////////////////////////////////////////////////////////////
  1123. BOOL CNLequation::GetRootsetNewton(int n, double x[], double t, double h, int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
  1124. { 
  1125.     int i,j,l;
  1126.     double am,z,beta,d,*y,*a,*b;
  1128.     y = new double[n];
  1130. // 构造矩阵
  1131. CMatrix mtxCoef(n, n);
  1132. CMatrix mtxConst(n, 1);
  1133.     a = mtxCoef.GetData();
  1134.     b = mtxConst.GetData();
  1136. // 迭代求解
  1137.     l=nMaxIt; 
  1138. am=1.0+eps;
  1139.     while (am>=eps)
  1140.     { 
  1141. Func(x,b);
  1142.         
  1143. am=0.0;
  1144.         for (i=0; i<=n-1; i++)
  1145.         { 
  1146. z=fabs(b[i]);
  1147.             if (z>am) 
  1148. am=z;
  1149.         }
  1150.         
  1151. if (am>=eps)
  1152.         { 
  1153. l=l-1;
  1154.             if (l==0)
  1155.             { 
  1156. delete[] y;
  1157.                 return FALSE;
  1158.             }
  1159.             
  1160. for (j=0; j<=n-1; j++)
  1161.             { 
  1162. z=x[j]; 
  1163. x[j]=x[j]+h;
  1164.                 
  1165. Func(x,y);
  1166.                 
  1167. for (i=0; i<=n-1; i++) 
  1168. a[i*n+j]=y[i];
  1169.                 
  1170. x[j]=z;
  1171.             }
  1172.             
  1173. // 调用全选主元高斯消元法
  1174. CLEquations leqs(mtxCoef, mtxConst);
  1175. CMatrix mtxResult;
  1176. if (! leqs.GetRootsetGauss(mtxResult))
  1177. {
  1178. delete[] y;
  1179. return FALSE;
  1180. }
  1182. mtxConst = mtxResult;
  1183. b = mtxConst.GetData();
  1185. beta=1.0;
  1186.             for (i=0; i<=n-1; i++) 
  1187. beta=beta-b[i];
  1188.             
  1189. if (beta == 0.0)
  1190.             { 
  1191. delete[] y;
  1192.                 return FALSE;
  1193.             }
  1194.             
  1195. d=h/beta;
  1196.             for (i=0; i<=n-1; i++) 
  1197. x[i]=x[i]-d*b[i];
  1198.             
  1199. h=t*h;
  1200.         }
  1201.     }
  1202.     
  1203. delete[] y;
  1204.     
  1205. // 是否在有效迭代次数内达到精度
  1206. return (nMaxIt>l);
  1207. }
  1209. //////////////////////////////////////////////////////////////////////
  1210. // 求非线性方程组最小二乘解的广义逆法
  1211. //
  1212. // 调用时，1. 须覆盖计算方程左端函数f(x)值及其偏导数值的虚函数
  1213. //             double Func(double x[], double y[])
  1214. //         2. 须覆盖计算雅可比矩阵函数的虚函数
  1215. //             double FuncMJ(int n, double x[], double y[])
  1216. //
  1217. // 参数：
  1218. // 1. int m - 方程的个数
  1219. // 2. int n - 未知数的个数
  1220. // 3. double x[] - 一维数组，长度为n，存放一组初值x0, x1, …, xn-1，要求
  1221. //                 不全为0，返回时存放方程组的最小二乘解，当m=n时，即是
  1222. //                 非线性方程组的解
  1223. // 4. double eps1 - 最小二乘解的精度控制精度，默认值为0.000001
  1224. // 5. double eps2 - 奇异值分解的精度控制精度，默认值为0.000001
  1225. //
  1226. // 返回值：BOOL 型，求解是否成功
  1227. //////////////////////////////////////////////////////////////////////
  1228. BOOL CNLequation::GetRootsetGinv(int m, int n, double x[], double eps1 /*= 0.000001*/, double eps2 /*= 0.000001*/)
  1229. { 
  1230. int i,j,k,l,kk,jt;
  1231.     double y[10],b[10],alpha,z,h2,y1,y2,y3,y0,h1;
  1232.     double *p,*d,*dx,*w;
  1233.     
  1234. // 控制参数
  1235. int ka = max(m, n)+1;
  1236.     w = new double[ka];
  1237.     
  1238. // 设定迭代次数为60，迭代求解
  1239. l=60; 
  1240. alpha=1.0;
  1241.     while (l>0)
  1242.     { 
  1243. CMatrix mtxP(m, n), mtxD(m, 1);
  1244. p = mtxP.GetData();
  1245. d = mtxD.GetData();
  1247. Func(x,d);
  1248.         FuncMJ(n,x,p);
  1250. // 构造线性方程组
  1251. CLEquations leqs(mtxP, mtxD);
  1252. // 临时矩阵
  1253. CMatrix mtxAP, mtxU, mtxV;
  1254. // 解矩阵
  1255. CMatrix mtxDX;
  1256. // 基于广义逆的最小二乘解
  1257. if (! leqs.GetRootsetGinv(mtxDX, mtxAP, mtxU, mtxV, eps2))
  1258.         { 
  1259. delete[] w;
  1260. return FALSE;
  1261.         }
  1262.         
  1263. dx = mtxDX.GetData();
  1265. j=0; 
  1266. jt=1; 
  1267. h2=0.0;
  1268.         while (jt==1)
  1269.         { 
  1270. jt=0;
  1271.             if (j<=2) 
  1272. z=alpha+0.01*j;
  1273.             else 
  1274. z=h2;
  1275.             
  1276. for (i=0; i<=n-1; i++) 
  1277. w[i]=x[i]-z*dx[i];
  1278.             
  1279. Func(w,d);
  1280.             
  1281. y1=0.0;
  1282.             for (i=0; i<=m-1; i++) 
  1283. y1=y1+d[i]*d[i];
  1284.             for (i=0; i<=n-1; i++)
  1285. w[i]=x[i]-(z+0.00001)*dx[i];
  1286.             
  1287. Func(w,d);
  1288.             
  1289. y2=0.0;
  1290.             for (i=0; i<=m-1; i++) 
  1291. y2=y2+d[i]*d[i];
  1292.             
  1293. y0=(y2-y1)/0.00001;
  1294.             
  1295. if (fabs(y0)>1.0e-10)
  1296.             { 
  1297. h1=y0; h2=z;
  1298.                 if (j==0) 
  1299. { 
  1300. y[0]=h1; 
  1301. b[0]=h2;
  1302. }
  1303.                 else
  1304.                 { 
  1305. y[j]=h1; 
  1306. kk=0; 
  1307. k=0;
  1308.                     while ((kk==0)&&(k<=j-1))
  1309.                     { 
  1310. y3=h2-b[k];
  1311.                         if (fabs(y3)+1.0==1.0) 
  1312. kk=1;
  1313.                         else 
  1314. h2=(h1-y[k])/y3;
  1315.                         
  1316. k=k+1;
  1317.                     }
  1318.                     
  1319. b[j]=h2;
  1320.                     if (kk!=0) 
  1321. b[j]=1.0e+35;
  1322.                     
  1323. h2=0.0;
  1324.                     for (k=j-1; k>=0; k--)
  1325. h2=-y[k]/(b[k+1]+h2);
  1326.                     
  1327. h2=h2+b[0];
  1328.                 }
  1329.                 
  1330. j=j+1;
  1331.                 if (j<=7) 
  1332. jt=1;
  1333.                 else 
  1334. z=h2;
  1335.             }
  1336.         }
  1337.         
  1338. alpha=z; 
  1339. y1=0.0; 
  1340. y2=0.0;
  1341.         for (i=0; i<=n-1; i++)
  1342.         { 
  1343. dx[i]=-alpha*dx[i];
  1344.             x[i]=x[i]+dx[i];
  1345.             y1=y1+fabs(dx[i]);
  1346.             y2=y2+fabs(x[i]);
  1347.         }
  1348.         
  1349. // 求解成功
  1350. if (y1<eps1*y2)
  1351.         { 
  1352. delete[] w;
  1353. return TRUE;
  1354.         }
  1355.         
  1356. l=l-1;
  1357.     }
  1358.     
  1359. delete[] w;
  1361. // 求解失败
  1362. return FALSE;
  1363. }
  1365. //////////////////////////////////////////////////////////////////////
  1366. // 求非线性方程组一组实根的蒙特卡洛法
  1367. //
  1368. // 调用时，须覆盖计算方程左端模函数值||F||的虚函数
  1369. //         double Func(int n, double x[])
  1370. //         其返回值为Sqr(f1*f1 + f2*f2 + … + fn*fn)
  1371. //
  1372. // 参数：
  1373. // 1. double x[] - 一维数组，长度为n，存放一组初值，返回时存放方程的一组实根
  1374. // 2. double xStart - 均匀分布的端点初值
  1375. // 3. int nControlB - 控制参数
  1376. // 4. double eps - 控制精度，默认值为0.000001
  1377. //
  1378. // 返回值：无
  1379. //////////////////////////////////////////////////////////////////////
  1380. void CNLequation::GetRootsetMonteCarlo(int n, double x[], double xStart, int nControlB, double eps /*= 0.000001*/)
  1381. { 
  1382. int k,i;
  1383.     double a,r,*y,z,z1;
  1385.     y = new double[n];
  1386.     
  1387. // 初值
  1388. a=xStart; 
  1389. k=1; 
  1390. r=1.0; 
  1392. z = Func(n, x);
  1394. // 用精度控制迭代求解
  1395.     while (a>=eps)
  1396.     { 
  1397. for (i=0; i<=n-1; i++)
  1398. y[i]=-a+2.0*a*rnd(&r)+x[i];
  1399.         
  1400. z1=Func(n,y);
  1401.         
  1402. k=k+1;
  1403.         if (z1>=z)
  1404.         { 
  1405. if (k>nControlB) 
  1406. { 
  1407. k=1; 
  1408. a=a/2.0; 
  1409. }
  1410. }
  1411.         else
  1412.         { 
  1413. k=1; 
  1414.             for (i=0; i<=n-1; i++) 
  1415. x[i]=y[i];
  1416.             
  1417. // 求解成功
  1418. z=z1;
  1419.             if (z<eps)  
  1420. {
  1421. delete[] y;
  1422. return;
  1423. }
  1424.         }
  1425.     }
  1427. delete[] y;
  1428. }