[原创]：用Brent方法和弦截法求一元方程的解！

highflybird

发一个一元方程的求解lisp程序。采用弦截法，对大多数方程有效。
附件中包含一个lisp程序和一个对话框，把对话框文件拷贝到CAD支持文件目录中（例如C:\Program Files\AutoCAD 2004\Support目录中）或者在支持文件搜索路径中添加指向这个对话框的目录，然后加载程序，运行solve,在对话框中输入数据运行即可。
[php]
(defun C:solve (/ f x0 n y id err result)
  (setq id (load_dialog "solve.dcl"))
  (if (new_dialog "dcl_solve" id)
    (progn
      (action_tile "p1" "(setq f $value)")
      (action_tile "p2" "(setq x0 $value)")
      (action_tile "p3" "(setq n $value)")
      (action_tile "accept" "(done_dialog)")
      (action_tile "help" "(choose 1)")
      (start_dialog)
      (unload_dialog id)
      (if (and x0 n f)
        (setq x0 (read x0) n (read n))
      )
      (if (or (= (type x0) 'REAL) (= (type x0) 'INT))
        (if (or (= (type x0) 'REAL) (= (type x0) 'INT))
          (progn
            (setq x0 (float x0))
            (setq n (fix (abs n)))
            (setq err (vl-catch-all-apply 'solve (list f x0 n)))
            (if (vl-catch-all-error-p err)
              (alert (strcat "方程无解或者" (vl-catch-all-error-message err)))
              (progn
                (alert result)
                (princ result)
              )
            )
          )
        )
      )
    )
  )
  (princ)
)
(defun diff1 (/ y0 y1 dy y2 x2)
  (setq y0 (fun f x0))                           ;第1次的函数值
  (setq y1 (fun f x1))                           ;第2次的函数值地
  (setq dy (- y1 y0))                            ;函数的差值
  (setq y2 (/ (- x1 x0) (- y1 y0)))              ;差分
  (setq x2 (- x1 (* y1 y2)))                     ;迭代得到新x值
  (setq x0 x1)                                   ;迭代第1次的x值
  (setq x1 x2)                                   ;迭代第2次的x值
)
(defun solve (f x0 n / dx x0 x1 x y i fra_x int_x)
  (VL-LOAD-COM)
  (if (> n 15)(setq n 15))                       ;如果超出精度就取15
  (setq dx 1e-2)                                 ;开始时的增量
  (setq x x0)                                    ;第1次的x值
  (setq x1 (+ x0 dx))                            ;第2次的x值
  (setq y (diff1))                               ;至少求1次迭代的值
  (setq i 1)                                     ;计数器，用来计算循环次数
  (while (and (> (abs (- y x)) (expt 0.1 (1+ n)));如果不满足精度
              (< i 1e9))                         ;且循环次数没有达到极限
    (setq x y y (diff1))                         ;则迭代
    (setq i (1+ i))                              ;计数器加1
  )
  (setq fra_x (rtos (- (abs x) (fix (abs x))) 2 n))
  (setq fra_x (vl-string-left-trim "0" fra_x))
  (setq int_x (rtos (fix x)))
  (setq result (strcat "\n方程" f "=0" "的解为:" int_x fra_x))
  ;;(setq result (strcat "\n方程" f "=0" "的解为:" (rtos y 2 n)))
)
(defun choose (n)
  (if (= n 1)
    (set_tile "error" "只接受x(小写)为变量,无须写等式,例如求解x^2-2=0写作x^2-2")
  )
)
;;;用方式1定义表达式求值函数
(defun fun (f x / f)
  (while (Wcmatch f "*x*")
    (setq f (vl-string-subst (rtos x 2 20) "x" f))
  )
  (vla-eval
    (vlax-get-acad-object)
    (strcat "ThisDrawing.SetVariable \"users1\",Cstr(" f ")")
  )
  (atof (getvar "users1"))
)
;;;用方式2定义表达式求值函数
(defun fun2 (f x) (arxload "geomcal.arx") (cal f))

[/php]

对话框文件：
[php]
//solve.dcl
dcl_solve : dialog {
  label = "解一元方程";
  key = "p0";
  : boxed_column {
    : row {
      key = "expression";
      : edit_box {
        key="p1";
        label= "请输入表达式:";
      }
    }
    : row {
      key = "F&A";
      fixed_width = true;
      : edit_box {
        key="p2";
        label= "初始值:";
        edit_width = 12;
      }
      : edit_box {
        key="p3";
        label= "精确到小数后位数:";
        edit_width = 4;
      }
    }
    spacer_1;
  }
  : row {
    fixed_width = true;
    alignment = centered;
    : ok_button : retirement_button {
      label = "确定";
      key = "accept";
      is_default = true;
    }
    spacer;
    spacer;
    : cancel_button : retirement_button {
      label = "取消";
      key = "cancel";
      is_cancel = true;
    }
    spacer;
    spacer;
    : help_button : retirement_button {
      label = "帮助(&H)";
      key = "help";
      is_help = true;
    }
  }
  errtile;
}  
[/php]
如果不能下载，请到下面网址:
http://www.mjtd.com/BBS/dispbbs. ... ID=58188&page=1
或者拷贝上面的代码段。
以前曾经用导数法（牛顿法）得出一个lisp程序，但出错较多，所以没发上来。
再过几天，发一个用Van Wijingaarden-Dekker-Brent方法求其区间解的lisp程序。
希望大家多多提出意见。

今天终于把那个Van Wijingaarden-Dekker-Brent方法的求解程序弄好了。
因为这个方法是C语言的，我是一字一句地把C语言翻译成lisp语言，很吃力，看来真的要多学点C语言了。
这个方法比弦截法有效得多，而且不易出错，推荐用这个方法。
lisp程序
[php]
;;;主程序
(defun C:solve (/ f x1 x2 n id err result)
  (setq id (load_dialog "solve1.dcl"))
  (if (new_dialog "dcl_solve" id)
    (progn
      (action_tile "p1" "(setq f $value)")  ;从对话框中得到表达式
      (action_tile "p2" "(setq x1 $value)") ;从对话框中得到下届
      (action_tile "p3" "(setq x2 $value)") ;从对话框中得到上届
      (action_tile "p4" "(setq n $value)")  ;从对话框中得到精度
      (action_tile "accept" "(done_dialog)");确定
      (action_tile "help" "(choose 1)")     ;帮助
      (start_dialog)
      (unload_dialog id)
      (if (and x1 x2 n f)
        (setq x1 (read x1) x2 (read x2) n (read n)))
      (if (or (= (type x1) 'REAL) (= (type x1) 'INT))
        (if (or (= (type x2) 'REAL) (= (type x2) 'INT))
          (if (or (= (type n) 'REAL) (= (type n) 'INT))
            (progn
              (setq n (fix (abs n)) x1 (float x1) x2 (float x2))
              (setq err (vl-catch-all-apply 'zbrent (list f x1 x2 n)))
              (if (vl-catch-all-error-p err)
                (alert (strcat "方程无解或者" (vl-catch-all-error-message err)))
                (progn
                  (alert result)
                  (princ result)
                )
              )
            )
          )
        )
      )
    )
  )
  (princ);很好看，我甚至不想用这个princ函数了:-)
)
;;;*********************************************
;;;用Van Wijingaarden-Dekker-Brent方法求解的函数
;;;*********************************************
(defun zbrent (f x1 x2 n / a b c d e fa fb fc EPS ITMAX i xm
                             tol tol1 min1 min2 p q r s fra_x int_x)
  (setq tol (expt 0.1 n))
  (if (> x1 x2)
    (setq a x2 b x1 c x1)
    (setq a x1 b x2 c x2)
  )
  (setq fa (func a))
  (setq fb (func b))
  (setq EPS 1e-16)
  (setq ITMAX 1e6)
  (if (or (and (> fa 0) (> fb 0)) (and (< fa 0) (< fb 0)))
    (princ "\nRoot must be bracketed in zbrent!")
  )
  (setq fc fb)
  (setq i 1)
  (while (<= i ITMAX)
  ;;(repeat 100
    (if (or (and (> fb 0) (> fc 0))
            (and (< fb 0) (< fc 0)))
      (setq c a fc fa d (- b a) e d);对a,b,c更名并调整解区间d
    )
    (if (< (abs fc) (abs fb))
      (setq a b b c c a fa fb fb fc fc fa)                     
    )
    (setq tol1 (+ (* 2.0 EPS (abs b)) (* 0.5 tol)));收敛性检查
    (setq xm (* 0.5 (- c b)))
    (if (or (<= (abs xm) tol1) (equal fb 0 1e-16));???????
      (setq i ITMAX)
    )
    (if (and (>= (abs e) tol1)
             (> (abs fa) (abs fb)))
      (progn
        (setq s (/ fb fa));将进行第二次反插
        (if (equal a c 1e-16)
          (setq p (* 2.0 xm s)
                q (- 1.0 s)
          )
          (setq q (/ fa fc)
                r (/ fb fc)
                p (* s (- (* 2.0 xm q (- q r)) (* (- b a) (1- r))))
                q (* (1- q) (1- r) (1- s))
          )
        )
        (if (> p 0) (setq q (- q)));检查是否在解区间内
        (setq p (abs p))
        (setq min1 (- (* 3.0 xm q) (abs (* tol1 q))))
        (setq min2 (abs (* e q)))
        (if (< (* 2.0 p) (if (< min1 min2) min1 min2))
          (setq e d d (/ p q));可以进行内插
          (setq d xm e d);不能进行内插，改用二分法
        )
      )
      (setq d xm e d);届下降速度太慢，改用二分法
    )
    (setq a b);将最新估值赋给a
    (setq fa fb)
    (if (> (abs d) tol1);计算新的实验解
      (setq b (+ b d))
      (setq b (+ b (sign tol1 xm)))
    )
    (setq fb (func b))
    (setq i (1+ i))
  )
  (if (equal (func b) 0 0.1)
    (setq fra_x (rtos (- (abs b) (fix (abs b))) 2 n)
        fra_x (vl-string-left-trim "0" fra_x)
        int_x (rtos (fix b))
        result (strcat "\n方程" f "=0" "的解为:" int_x fra_x)
    )
    (setq result (strcat "\n方程" f "=0" "在此区间无解"))
  )
)
(defun sign (x y)
  (if (>= y 0)
    (abs x)
    (- (abs x))
  )
)
(defun choose (n)
  (if (= n 1)
    (set_tile "error" "    方程只接受x(小写)为变量,无须写等式,例如求解x^2-2=0写作x^2-2")
  )
)
;;;用方式1定义表达式求值函数
(defun func1 (x)
  (arxload "geomcal.arx")
  (cal f)
)
;;;用方式2定义表达式求值函数
(defun func (x)
  (setq fun f)
  (while (Wcmatch fun "*x*")
    (setq fun (vl-string-subst (rtos x 2 20) "x" fun))
  )
  (vla-eval
    (vlax-get-acad-object)
    (strcat "ThisDrawing.SetVariable \"users1\",Cstr(" fun ")")
  )
  (read (getvar "users1"))
)
[/php]
对话框文件
[php]
//solve.dcl
dcl_solve : dialog {
  label = "一元方程求解程序";
  key = "p0";
  : boxed_column {
    : row {
      key = "expression";
      : edit_box {
        key="p1";
        label= "一元方程:";
      }
    }
    : row {
      key = "F&A";
      fixed_width = true;
      : edit_box {
        key="p2";
        label= "区间下届:";
        edit_width = 8;
      }
      : edit_box {
        key="p3";
        label= "区间上届:";
        edit_width = 8;
      }      
      : edit_box {
        key="p4";
        label= "精确位数:";
        edit_width = 2;
      }
    }
    spacer_1;
  }
  ok_cancel_help;
  errtile;
}  
[/php]
下面是效果截图：

                               
登录/注册后可看大图

用计算器验证了程序得出的值，完全满足要求。

highflybird

稍微做了不得一下修改，修正了几个bug,
重新发上来:
[php]

;;;主程序
(prompt "请输入solve命令!")
(defun C:solve (/ f x1 x2 n id err result str1 str2 str3 str4 str5)
  (setq id (load_dialog "solve1.dcl"))
  (if (new_dialog "dcl_solve" id)
    (progn
      (action_tile "p1" "(setq f $value)")  ;从对话框中得到表达式
      (action_tile "p2" "(setq x1 $value)") ;从对话框中得到下届
      (action_tile "p3" "(setq x2 $value)") ;从对话框中得到上届
      (action_tile "p4" "(setq n $value)")  ;从对话框中得到精度
      (action_tile "help" "(choose 1)")     ;帮助
      (start_dialog)
      (unload_dialog id)
      (if (and x1 x2 n f)
        (progn
          (setq x1 (atof x1) x2 (atof x2) n (abs (atoi n)) n (if (> n 20) 20 n))
          (setq err (vl-catch-all-apply 'zbrent (list f x1 x2 n)))
          (if (vl-catch-all-error-p err)
            (alert (strcat "方程无解或者" (vl-catch-all-error-message err)))
            (progn
              (setq str1 (car result)
                    str2 "\n下面是求解验证结果:"
                    str3 (cadr result)
                    str4 (caddr result)
                    str5 (strcat str1 str2 "\nf(" str3 ")=" str4)
              )
              (alert str5)
              (princ str5)
            )
          )
        )
      )
    )
  )
  (princ)
)
;;;*********************************************
;;;用Van Wijingaarden-Dekker-Brent方法求解的函数
;;;*********************************************
(defun zbrent (f x1 x2 n / a b c d e fa fb fc EPS ITMAX i xm test
                           tol tol1 min1 min2 p q r s fra_x int_x)
  (setq tol (expt 0.1 n));容差
  (setq ITMAX 10000)     ;最大迭代次数
  (setq EPS 1e-16)       ;计算机有效精度
  (if (> x1 x2)          ;如果下届大于上届
    (setq a x2 b x1 c x1);则交换上下界
    (setq a x1 b x2 c x2)
  )
  (setq fa (func a))
  (setq fb (func b))
  (if (or (and (> fa 0) (> fb 0)) (and (< fa 0) (< fb 0)));上下界函数值同号
    (princ "\nRoot must be bracketed in zbrent!")
  )
  (setq fc fb)
  (setq i 1)
  (while (<= i ITMAX)
    (if (or (and (> fb 0) (> fc 0))
            (and (< fb 0) (< fc 0)))
      (setq c a fc fa d (- b a) e d);对a,b,c更名并调整解区间d
    )
    (if (< (abs fc) (abs fb))
      (setq a b b c c a fa fb fb fc fc fa)                     
    )
    (setq tol1 (+ (* 2.0 EPS (abs b)) (* 0.5 tol)));收敛性检查
    (setq xm (* 0.5 (- c b)))
    (if (or (<= (abs xm) tol1) (equal fb 0 1e-16));跳出循环
      (setq i ITMAX)
      (progn
        (if (and (>= (abs e) tol1)
                 (> (abs fa) (abs fb)))
          (progn
            (setq s (/ fb fa));将进行第二次反插
            (if (equal a c 1e-16)
              (setq p (* 2.0 xm s)
                    q (- 1.0 s)
              )
              (setq q (/ fa fc)
                    r (/ fb fc)
                    p (* s (- (* 2.0 xm q (- q r)) (* (- b a) (1- r))))
                    q (* (1- q) (1- r) (1- s))
              )
            )
            (if (> p 0) (setq q (- q)));检查是否在解区间内
            (setq p (abs p))
            (setq min1 (- (* 3.0 xm q) (abs (* tol1 q))))
            (setq min2 (abs (* e q)))
            (if (< (* 2.0 p) (if (< min1 min2) min1 min2))
              (setq e d d (/ p q));可以进行内插
              (setq d xm e d);不能进行内插，改用二分法
            )
          )
          (setq d xm e d);届下降速度太慢，改用二分法
        )
        (setq a b);将最新估值赋给a
        (setq fa fb)
        (if (> (abs d) tol1);计算新的实验解
          (setq b (+ b d))
          (setq b (+ b (if (>= xm 0)(abs tol1) (- (abs tol1)))))
        )
        (setq fb (func b))
      )
    )
    (setq i (1+ i))
  )
  (setq test (func b))
  (setq fra_x (rtos (- (abs b) (fix (abs b))) 2 n))
  (setq fra_x (vl-string-left-trim "0" fra_x))
  (setq int_x (itoa (fix b)))
  (if (equal test 0 1)
    (setq result (strcat "\n方程" f "=0" "的解为:" int_x fra_x)
          result (list result (strcat int_x fra_x) (rtos test 1 n))
    )
    (setq result (strcat "\n方程" f "=0" "在此区间无解")
          result (list result (strcat int_x fra_x) (rtos test 1 n))
    )
  )
)
;;;帮助说明函数
(defun choose (n)
  (if (= n 1)
    (alert "方程式只接受x(小写)为变量,不规范很可能出错!
           \n无须写等式,例如求解x^2-2=0写作x^2-2.
           \n程序采用brent方法求解,不保证每个方程都有效!
           \n有什么问题email: highflybird@sina.com"
    )
    ;;(set_tile "error" "方程式只接受x(小写)为变量,无须写等式,例如求解x^2-2=0写作x^2-2.")
  )
)
;;;用方式1定义表达式求值函数
(defun func1 (x)
  (arxload "geomcal.arx")
  (cal f)
)
;;;用方式2定义表达式求值函数
(defun func (x / fun vbeval)
  (setq fun f)
  (while (Wcmatch fun "*x*")
    (setq fun (vl-string-subst (rtos x 2 20) "x" fun))
  )
  (vla-eval
    (vlax-get-acad-object)
    (strcat "ThisDrawing.SetVariable \"users1\",Cstr(" fun ")")
  )
  (atof (getvar "users1"))
)

[/php]
对话框文件:
[php]
//solve.dcl
dcl_solve : dialog {
  label = "一元方程求解程序";
  key = "p0";
  : boxed_column {
    : row {
      key = "expression";
      : edit_box {
        key="p1";
        label= "一元方程:";
      }
    }
    : row {
      key = "F&A";
      fixed_width = true;
      : edit_box {
        key="p2";
        label= "区间下届:";
        edit_width = 8;
      }
      : edit_box {
        key="p3";
        label= "区间上届:";
        edit_width = 8;
      }      
      : edit_box {
        key="p4";
        label= "精确位数:";
        edit_width = 2;
      }
    }
    spacer_1;
  }
  ok_cancel_help;
  //errtile;
}  
[/php]