数值计算之解方程

此处用比较纯LISP的方法来解非线性方程。
主要采用Van Wijingaarden-Dekker-Brent方法求解。

下面是其源码:

;;;*********************************************
;;;用Van Wijingaarden-Dekker-Brent方法求解的函数
;;;*********************************************
(defun zbrent (f x1 x2 n / a b c d e i p q r s xm fa fb fc EPS Itr tol tol1 min1 min2)
  (setq tol (expt 0.1 n))					;容差
  (setq Itr 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 Itr)
    (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 Itr)
      (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))
  )
  b
)

下面的代码包含了一个应用及其对话框例子:

(vl-load-com)
;;;主程序
(prompt "请输入solve命令!")
 
;;;表达式求值
(defun Getfunc (str)
  (CAL:Expr2Func str 'func '(x))
)
 
(defun C:solve (/ f x1 x2 n id Dcl_File)
  (setq id (load_dialog (setq Dcl_File (Write_Dcl))))		;装入对话框
  (vl-file-delete Dcl_File)					;删除临时文件
  (if (new_dialog "dcl_solve" id)
    (progn
      (action_tile "EXP" "(Getfunc (setq f $value))") 		;从对话框中得到表达式
      (action_tile "MIN" "(setq x1 $value)") 			;从对话框中得到下届
      (action_tile "MAX" "(setq x2 $value)") 			;从对话框中得到上届
      (action_tile "PRC" "(setq n $value)") 			;从对话框中得到精度
      (action_tile "ANS" "(GetAnswer f x1 x2 n)")               ;求解
      (action_tile "help" "(choose 1)")				;帮助
      (start_dialog)
      (unload_dialog id)
      (princ (GetAnswer f x1 x2 n))
    )
  )
  (princ)
)
 
;;;*********************************************
;;;从对话框得到方程参数并求解和显示结果         
;;;*********************************************
(defun GetAnswer (f x1 x2 n / str1 str2 str3 str4 str5 ret test fra_x int_x result)
  (if (and x1 x2 n f)
    (progn
      (setq x1 (atof x1))
      (setq x2 (atof x2))
      (setq n  (abs (atoi n)))
      (if (> n 20)
         (setq n 20)
      )
      (GetFunc f)
      (setq ret (vl-catch-all-apply 'zbrent (list f x1 x2 n)))
      (if (vl-catch-all-error-p ret)
	(setq str5 (strcat "出错: " (vl-catch-all-error-message ret)))
	(progn
	  (setq test (func ret))
	  (setq fra_x (rtos (- (abs ret) (fix (abs ret))) 2 n))
	  (setq fra_x (vl-string-left-trim "0" fra_x))
	  (setq int_x (itoa (fix ret)))
	  (if (equal test 0 1)
	    (setq result (strcat "方程" f "=0" "的解为:" int_x fra_x)
		  result (list result (strcat int_x fra_x) (rtos test 2 n))
	    )
	    (setq result (strcat "方程" f "=0" "在此区间无解")
		  result (list result (strcat int_x fra_x) (rtos test 2 n))
	    )
	  )
	  (setq	str1 (car result)
		str2 "\n下面是求解验证结果:"
		str3 (cadr result)
		str4 (caddr result)
		str5 (strcat str1 str2 "\nf(" str3 ")=" str4)
	  )
	)
      )
      (set_tile "error" str5)
      str5
    )
    (set_tile "error" "输入有误!")
  )
)
 
 
 
;;;帮助说明函数
(defun choose (n)
  (if (= n 1)
    (alert
      "方程式只接受x(小写)为变量,不规范很可能出错!
           \n无须写等式,例如求解x^2-2=0写作x^2-2.
           \n程序采用brent方法求解,不保证每个方程都有效!
	   \n有什么问题email: highflybird@qq.com"
    )
    (set_tile
      "error"
      "方程式只接受x(小写)为变量,无须写等式,例如求解x^2-2=0写作x^2-2."
    )
  )
)
 
;;; for DCL
(defun Write_Dcl (/ Dcl_File file str)
  (setq Dcl_File (vl-filename-mktemp nil nil ".DCL"))
  (setq file (open Dcl_File "w"))
  (princ
    "dcl_solve : dialog {
      label = \"一元方程求解程序\";
      : boxed_column {
        : edit_box {
          key=\"EXP\";
          label= \"一元方程:\";
        }
        : row {
          : edit_box {
            key=\"MIN\";
	    label= \"区间下届:\";
            edit_width = 8;
          }
          : edit_box {
            key=\"MAX\";
            label= \"区间上届:\";
            edit_width = 8;
          }
          : edit_box {
            key=\"PRC\";
            label= \"精确位数:\";
            edit_width = 2;
          }
        }
        spacer_1;
      }
      : row {
        : button {
          key=\"ANS\";
          label=\"求解\";
        }
        ok_cancel_help;
      }
      errtile;
    }"
    file
  )
  (close file)
  Dcl_File
)

求解速度还可以,对于一般的,如果估算好求解区间,都会得到正确结果。


发表回复

您的电子邮箱地址不会被公开。 必填项已用 * 标注