此处用比较纯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
)求解速度还可以，对于一般的，如果估算好求解区间，都会得到正确结果。