记忆从这里开始

时间过去了很久,希望有些东西能沉淀下来,虽然我有记录的爱好,却总不能坚持。 陆陆续续这么多年,弄了很多很多博客,要不是被限制,要不就是被墙了,甚至无缘无故被删除掉,白白浪费了我很多心血。所以我决定在建军节前夕,购买了这个空间和域名,在此安家。 这是博客的第一篇,特此记录。

CAD问题小百科

我是搬运工:

从这里摘录了一些常见的CAD小问题:

https://www.cnblogs.com/JJBox/p/10848766.html

一、词典问题集合
这是一个常见的问题。
先说解决方案:
方法1:可以用删除词典的方法: 

(dictremove(namedobjdict)"ACAD_DGNLINESTYLECOMP")

方法2:
如果你有用源泉建筑插件可以用pua命令也可以用桌子带有的工具,
点我
方法3:
另存R14格式,再保存成你的版本。
关闭之后再次打开,保存,速度变快。(尝试清理词典)

高版本貌似pu就有选项了就可以了…

这个问题症状如下:
1、CAD保存慢:
2、cad很多线型!(尝试清理词典)CAD不能粘贴,出现_pasteclip:
3、dgn线型的问题,天正也会导致这个问题(尝试清理词典)
4、在无天正的CAD状态下打开含有天正图元的DWG,装上天正,然后转旧版本;
5、CAD每次保存都会出现这个: *警告* 多重从属对象,句柄“FXXXXX”;(尝试清理词典)

二、注册问题

cad 打开就闪退 用Everything搜索 删除 flexnet 目录下所有文件,这是由于注册文件导致的…需要重新注册…
cad2012~CAD2017打开慢,在注册表编辑器中找到下面位置:
HKEY_CURRENT_USER\Software\Autodesk\AutoCAD\R18.2\ACAD-A005:804\InfoCenter\InfoCenterOn 然后把值改为0
也可能是默认打印机设置为网络打印机,网络打印机关机了,检索不到,造成缓慢…代码的结局方案可以看我另一篇代码解决网络打印机的问题:博文

三、权限问题
win10禁止拖拉程序到绘图区?
运行regedit,注册表 HKEY_LOCAL_MACHINE\SOFTWARE\Microsoft\Windows\CurrentVersion\Policies\System\EnableLUA 
然后把值改为0
但是我告诉你…会导致win10的商店用不了..貌似内置应用也用不了..你可能是正版家庭版win用户,建议更换专业版。。。
win10的cad用ctrl+q会是_copyclip
把一个cad改成兼容性win7执行就可以了…但是貌似治标不治本,形而上学啊….万一兼容性会导致某些问题呢….
你可能是正版受害者…..
CAD搜索报错,高版本ctrl+u弹界面
CAD因为会搜索当前dwg路径的所有文件,如果dwg在桌面,那么快捷方式无效会导致报错.
说实话没找到解决方案,桌子为什么那么蛋疼做这个东西…

四、CAD界面跳动
Statusbar,autowrap,0 ;2018关闭命令栏跳动(跳动也可以去关闭坐标显示)

五、CAD操作问题
属性名称不能相同,做成属性块之后会红字显示,并且用ATTSYNC更新时候会出错,拉伸属性的时候会出错.
有时候在位编辑之后,保存在位编辑了,可是参照面板无法动了,这个时候只要保存一下,然后ctrl+z,就可以解决了..蛮神奇的…
代码的话看我另一篇win32打开参照面板的博文的主函数上面,
原理就是发送命令,无他…有知道原因的告诉我一下..

六、复制粘贴问题
症状1: cad2008复制粘贴就崩溃但是高版本不会,可能是天正电气的一个bug导致的.
症状2: 之前也试过,用高版本的CAD清理了之后,就可以,画了一些东西之后就又不行了.
解答: 要重新用最新版的天正和其匹配的cad进行转旧出来.

七、CAD安装问题
cad2007遇到了这种情况
安装问题安装CAD出现C++2005问题的解决方法,出现此问题,
原因大部分就是因为安装路径目录或是下载的安装包有中文字体,导致电脑不能识别,只要把中文该成英文字母形式就可以了。
详见: cad2012安装的问题
cad2020安装时候出现1603安装失败:
运行,cmd,输入:”C:\Program Files (x86)\Common Files\Autodesk Shared\AdskLicensing\Current\AdskLicensingService\AdskLicensingService.exe”
CAD2020激活界面问题 许可管理器不起作用或未正确安装
混合电源设置会影响Autodesk 2020产品许可
下载里面的Autodesk Licensing Service v9.0.3.46.zip 
安装即可.
20190713:这个方法并不是很好,似乎是我使用了最新的win1903的最新更新导致? 因为我今天打开又遇到了…
20191205:开启此服务,注意是不是你的电脑管家或者360优化关闭了相关服务.
右击我的电脑-管理-服务与应用程序-服务

删除所有注册表,再安装开始-运行-输入“regedit”,打开注册表,找到下面的注册表路径:HKEY_LOCAL_MACHINE\SOFTWARE\Autodesk\AutoCAD\R18.0
右击鼠标删除 R18.0 整个键值。
找到下面的注册表路径:
HKEY_CORRENT_user\SOFTWARE\Autodesk\AutoCAD\R18.0
右击鼠标删除 R18.0 整个键值。
最后一步,也是最重要的,找到下面的注册表路径:
HKEY_LOCAL_MACHINE\SOFTWARE\Classes\installer\Products 的下面,
右击鼠标删除“7D2F38751008某某某”开头的键值,删除以后就可以重新安装了 

八、图块编辑问题
很危险的问题****因为进块编辑器,编辑之后就退出,结果一卡,这个块把模型原有的图替代了,会丢失整个模型块表…
你的图要重画咯~~~~~~
画图记得备份,备份器c#的代码实现我公开在这里
某些人: Acad2014和浩辰2018…. 
瓦蓝:     win7下Acad2012字体问题集合

九、字体问题
cad 两台电脑的宋体相同,但是一台电脑保存另一台出现问号,删除KT_ST.ttf,
更为详尽的说明代码的解决方案
cad 高版本面板乱码,参照面板弹出之后,其他面板会乱码,删除所有cad支持目录的fonst文件夹中的.ttf  .ttc之类的字体,
包括插件和cad自带目录,高版本不知道为什么桌子自带这样的问题,真的令人头痛…..
插入或者绑定外部参照出现文字偏移,请不要使用”新宋体.TTF”,
详见点我 
什么都不加载,启动cad时也乱码: 系统问题,区域改成新加坡重启后,可以加载了,再改回中国,应当就ok了。

十、多窗口问题
CAD一个程序中打开多个窗口设置:
1、命令行输入SDI,回车,输入0,回车。
2、点击工具条上“窗口”,再点击下面的选项。
3、选择不同的在同一窗口打开多个文件的方式。
4、打开cad工具栏中:工具-选项-系统-(去掉单文档兼容模式前面的勾),然后就好了。
5、在工具栏的空白处点右键,然后在“选项卡”前面打上对勾,就可以在同一个cad程序里面打开多个cad文件了。

十一、打印问题
打印cad用pdf factory的时候发生打印错误,
提示:内部错误:PrepareDispatcher was unable to launch the
发生在联网打印就不行,断网打印就可以…官网的bat我也下载了,也不行…….
解决:我先是升级win10到1903版本…..
升级工具
开心了一个小时之后,我做了局域网共享的操作之后,发现又出现了!所以我肯定了是共享文件夹的方式出现了问题,
经过查询之后,我发现了要设置防火墙允许通过的列表中:远程桌面,远程协助,远程关机,
这几个要关掉.貌似就是这几个导致的(只是貌似,因为我先运行了下面的.bat)文件和打印机共享、网络发现、SMBDirect上的文件和打印机共享这几个要开启.
之前这个仁兄写的是关闭防火墙的,我建议他打开防火墙,所以你们可以用这个工具的最新版来解决局域网共享:
局域网共享一键修复.bat

十二、杂项
1、cad 菜单选择打开文件的时候,等半天才出现文件打开对话框,点我
2、cad 复制中断产生的不可见图元,无法刷新,点我
3、怪事1,CAD在局域网联机工作,一直可以保存的,但是突然变成只读,但是直接关掉CAD会丢失了文件!
怪事2:CAD在本机工作或联机工作,保存出现错误,生成了文件,但是直接关掉CAD会丢失了文件!
以上似乎是win7的问题,局域网共享的问题,详细我已经没有测试条件了…建议换win10
在企业版win中,我发现了加载了自己的插件,然后闪退,致命错误…这就很不解了….换专业版
在家庭版win中,使用源泉,貌似会发生不可思议的加载失效事情….换专业版
4、选项卡预览关闭
将FILETABPREVIEW系统变量的值设为0,将显示方式变为列表,无法关闭预览
5、autocad打开时多个文件时任务栏也同时出现多个窗口的解决办法:
打开CAD; 输入taskbar命令;CAD2011默认为1,改为0即可。

小问题,大思考

一个博士朋友的儿子刚上初中不久,他就提出了下面的一个问题,我感觉很有深度,很难得,现在贴出来说说。

如下面图所示:已知线段AB和一点P,n等分线段,得到n+1个点与P连线。

问:n趋于无穷大时,这些连线长的平均值趋向一个定数吗? 如果是,该怎样求?

粗看起来好像没什么难度,然而稍微计算一下,就发现不那么简单,博士拿出了Maple这个数学软件,用求极限法得到了一个并不简单的结果:

后来,我觉得这个极限可以化为积分,得到了相同的结果:

其中的字母含义如下:

这个结果显然不是初中的知识了。

这个小孩很有天赋,他还把这个问题做引申发散开来。

首先,他问这样的平均距离除了对直线有外,对圆也有么?

我开始用积分算了一下,得到一个巨复杂的结果:

后来稍微优化了一下,

但可以看出,这已经不是可以用初等函数来表达的了,已经用了椭圆积分来表达了。然后,他又继续问了一个问题,在一个三角形内,是否存在这样的一个点,使得点到这三条线段的上述平均值相等?

在我看来,一般情况下存在这个点,原因如下:
因为一点(x,y)到线段的平均值的积分解是一个与x,y有关的二元函数
因此 另   p1p2的 值  设为f1(x,y), p2p3的值设置为f2(x,y),p3p1的值设置为f3(x,y)
显然有  f1(x,y)=f2(x,y); f2(x,y)=f3(x,y)
这是一个二元方程组,应该有一个解,但 估计无法用公式求解,只能数值求解了。

下面贴出一些LISP代码来验证求解。

(defun ALG:AverageLength (p p1 p2 / A A1 A2 A3 A4 H HH L RT V)
  (setq L (distance p1 p2))
  (setq v (trans (mapcar '- p p1) 0 (mapcar '- p2 p1)))
  (setq h (abs (car v)))
  (setq a (- (caddr v)))
  (if (> L 0)
    (progn
      (setq a (/ a L))
      (setq h (/ h L))
      (setq hh (* h h))
      (setq a1 (+ (* a a) hh))
      (setq a2 (+ a1 a a 1))
      (setq a3 (sqrt a1))
      (setq a4 (sqrt a2))
      (if (equal hh 0 1e-8)
	(setq rt (+ (* a (- a4 a3)) a4))
        (setq rt (+ (* a (- a4 a3))
		    (* hh (- (log (+ a 1 a4)) (log (+ a a3))))
		    a4
	         )
        )
      )
      (setq rt (* 0.5 rt L))
    )
    (distance p p1)
  )
)
 
;(1/2)*a*sqrt(a^2+h^2)-(1/2)*h^2*ln(-a+sqrt(a^2+h^2))+(1/2)*b*sqrt(b^2+h^2)+(1/2)*h^2*ln(b+sqrt(b^2+h^2))
;;;采用积分法
(defun ALG:AverageLength2 (p p1 p2 / A AH B BH H HH L V rt)
  (setq L (distance p1 p2))
  (setq v (trans (mapcar '- p p1) 0 (mapcar '- p2 p1)))
  (setq h (abs (car v)))
  (setq a (caddr v))
  (setq b (- L a))
  (if (equal L 0 1e-8)
    (distance p p1)
    (if (equal h 0 1e-8)
      (/ (+ (* a a) (* b b)) 2 L)
      (setq a  (/ a L)
	    b  (/ b L)
	    h  (/ h L)
	    hh (* h h)
	    ah (sqrt (+ (* a a) hh))
	    bh (sqrt (+ (* b b) hh))
	    rt (* hh (log (/ (+ b bh) (- ah a))))
            rt (* 0.5 L (+ RT (* a ah) (* b bh)))
      )
    )
  )
)
 
;;;收敛很慢
(defun ALG:AverageLength1 (p p1 p2 / A AV D H HH I L N S UN V X)
  (setq L (distance p1 p2))
  (setq v (trans (mapcar '- p p1) 0 (mapcar '- p2 p1)))
  (setq h (abs (car v)))
  (setq a (- (caddr v)))
  (if (> L 0)
    (progn
      (setq n 1000)
      (setq i 0)
      (setq s 0)
      (setq hh (* h h))
      (setq un (/ L n))
      (repeat (1+ n)
	(setq x (+ a (* i un)))
        (setq d (sqrt (+ (* x x) hh)))
	(setq s (+ s d))
	(setq i (1+ i))
      )
      (setq av (/ s (1+ n)))
    )
    (distance p p1)
  )
)
 
(defun ALG:AverageLengthOfCircle (p c r / A D I L N U X Y)
  (setq p (mapcar '- p c))
  (setq i 0)
  (setq l 0)
  (setq n 36000)
  (setq u (/ pi n))
  (repeat n
    (setq a (/ i u))
    (setq x (* r (cos a)))
    (setq y (* r (sin a)))
    (setq d (distance p (list x y)))
    (setq l (+ d l))
    (setq i (1+ i))
  )
  (/ l n)
)
 
(defun c:xx ()
  (initget 9)
  (setq p (getpoint "\n直线外一点:"))
;;;  (initget 9)
;;;  (setq c (getpoint "\n圆心:"))
;;;  (initget 15)
;;;  (setq r (getdist C "\n半径:"))
  (setq c '(0 0))
  (setq r 1)
  (setq z (ALG:AverageLengthOfCircle p c r))
  (princ "\n所求的值为: ")
  (princ (rtos z 2 20))
  (princ)
)
 
(defun c:gg ()
  (setq i 0)
  (setq p0 '(0 0))
  (setq l 0)
  (setq x 0)
  (setq n 1000000)
  (setq u (/ 1.0 n))
  (repeat n
    (setq y (cos x))
    (setq d (distance p0 (list x y)))
    (setq l (+ l d))
    (setq x (* i u))
    (setq i (1+ i))
  )
  (setq v (/ l n))
  (princ "\n所求的值为: ")
  (princ (rtos v 2 20))
  (princ)
)
 
(defun c:ccc()
  (initget 9)
  (setq p (getpoint "\n直线外一点:"))
  (setq p (trans p 1 0))
  (initget 9)
  (setq p1 (getpoint "\n直线第一点:"))
  (setq p1 (trans p1 1 0))
  (initget 9)
  (setq p2 (getpoint "\n直线第一点:"))
  (setq p2 (trans p2 1 0))
  (setq x (ALG:AverageLength p p1 p2))
  (princ (strcat "\n所求的平均值是:" (rtos x 2 20)))
  (setq x1 (ALG:AverageLength1 p p1 p2))
  (princ (strcat "\n极限法所求的平均值是:" (rtos x1 2 20)))
 
  (setq x2 (ALG:AverageLength2 p p1 p2))
  (princ (strcat "\n积分法所求的平均值是:" (rtos x2 2 20)))
 
;;;  (UTI:BENCH 1000
;;;    (list
;;;      (list 'ALG:AverageLength p p1 p2)
;;;      (list 'ALG:AverageLength1 p p1 p2)
;;;    )
;;;  )
 
  (princ)
)

有兴趣的不妨测试一下。

数值计算之解方程

此处用比较纯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
)

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

数值计算之求积分

用LISP编写了一个求积分的程序:

里面采用了各种方法求积分和各种类型的积分。下面我把各种方法的源码贴出。

方法一: 勒让德-高斯积分法。

;;;=============================================================
;;; 计算勒让德-高斯求积函数的系数 ,如下面的6次项系数           
;;; p0(x) = 1                                                   
;;; p1(x) = x                                                   
;;; p2(x) = (3*x^2-1)/2                                         
;;; p3(x) = (5*x^3-3*x)/2                                       
;;; p4(x) = (35*x^4-30*x^2+3)/8                                 
;;; p5(x) = (63*x^5-70*x^3+15*x)/8                              
;;; p6(x) = (231*x^6-315*x^4+105*x^2-5)/16                      
;;; x-2                                                         
;;; x-1                                                         
;;; 9x-5                                                        
;;; 216*x^2-216*x+49                                            
;;; 45000*x^2-32200*x+5103                                      
;;; 2025000*x^3-2025000*x^2+629325*x-58564                      
;;; 142943535000*x^3-1130712534000*x^2+27510743799*x-1976763932 
;;; 输入: x1,x2 区间(一般来说是-1..1),n 迭代次数,eps迭代精度   
;;; 输出: 勒让德-高斯求积函数的系数,用点表集表示                
;;; http://mathworld.wolfram.com/Legendre-GaussQuadrature.html  
;;;=============================================================
(defun Math::Int:Legendre_Polynomial (x1 x2 n eps / xi wi FI I ITER J M P1 P2 P3 PP XL XM Z Z1)
  (setq m (/ (1+ n) 2))
  (setq xm (* 0.5 (+ x2 x1)))
  (setq xl (* 0.5 (- x2 x1)))
  (setq i 1)
  (repeat m
    (setq z (cos (/ (* pi (- i 0.25)) (+ n 0.5))))
    (setq iter 0)
    (while (and (not (equal z z1 eps)) (< iter 1000))
      (setq p1 1.0)
      (setq p2 0.0)
      (setq j 1)
      (repeat n
	(setq p3 p2)
	(setq p2 p1)
	(setq p1 (/ (- (* z p2 (+ j j -1.)) (* (1- j) p3)) j))
	(setq j (1+ j))
      )
      (setq pp (/ (* n (- (* z p1) p2)) (1- (* z z))))
      (setq z1 z)
      (setq z (- z1 (/ p1 pp)))
      (setq iter (1+ iter))
    )
    (setq fi (/ xl 0.5 (- 1 (* z z)) pp pp))
    (setq xi (cons (cons (1- i) (- xm (* xl z))) xi))
    (setq wi (cons (cons (1- i) fi) wi))
    (if (/= (1- i) (- n i))
      (setq xi (cons (cons (- n i) (+ xm (* xl z))) xi)
	    wi (cons (cons (- n i) fi) wi)
      )
    )
    (setq i (1+ i)) 
  )
  (MATH::INT:Bind xi wi)
)
 
;;;=============================================================
;;; 勒让德-高斯求积函数                                         
;;;=============================================================
(defun MATH::INT:Gauss-Legendre (a b eps / AA BB EP FX G H I L M P S W X)
  (setq l '((-0.93246951420315202787 . 0.17132449237917034504)  
	    (-0.66120938646626451363 . 0.36076157304813860756)     
	    (-0.23861918608319690859 . 0.46791393457269104739)  
	    ( 0.23861918608319690859 . 0.46791393457269104739)     
            ( 0.66120938646626451363 . 0.36076157304813860756)     
	    ( 0.93246951420315202787 . 0.17132449237917034504)
	  )
  )                                                         
  ;(setq L (Math::Int:Legendre_Polynomial -1 1 100 2e-20))
  (setq m 1)
  (setq h (- b a))
  (setq s (abs (* 0.001 h)))
  (setq p 1e100) 	
  (setq ep (1+ eps))
  (while (and (>= ep eps) (> (abs h) s))
    (setq g 0)
    (setq i 1)
    (repeat m
      (setq bb (+ a (* i h)))
      (setq aa (- bb h))
      (setq w 0)
      (foreach k l
	(setq x (* 0.5 (+ bb aa (* (- bb aa) (car k)))))
	(setq fx (MATH::INT:func x))
	(setq w (+ w (* fx (cdr k))))
      )
      (setq g (+ g w))
      (setq i (1+ i))
    ) 
    (setq g (* g h 0.5))
    (setq ep (/ (abs (- g p)) (1+ (abs g))))
    (setq p g)
    (setq m (1+ m))
    (setq h (/ (- b a) m))
  )
  g
)
 
;;;=============================================================
;;; 勒让德-高斯求积函数(另一方法,慢些)                         
;;;=============================================================
(defun MATH::INT:Gauss-Legendre1 (a b eps / FLAG G H L N X Y)
  (setq n 1)                                                      
  (setq flag T)							;是否进行迭代
  (while (and (< n 100) flag)	
    (setq g 0)
    (setq L (Math::Int:Legendre_Polynomial a b n eps))
    (foreach w L
      (setq x (car w))
      (setq y (MATH::INT:func x))
      (setq g (+ g (* (cdr w) y)))
    )
    (if (equal g h eps)
      (setq flag nil)
      (setq n (+ n n)
	    h g
      )
    )
  )
  g
)

方法二:高斯-埃米尔特积分

;;;=============================================================
;;; 高斯-埃米尔特积分                                           
;;; 功能: 计算 e^(-x^2)*f(x)的广义积分(在区间-INF..INF上)       
;;;=============================================================
(defun Math::INT:Gauss-Hermite (a b eps / n L g x y)
  (setq n 100)
  (setq L (Math::INT:GetHermite n))
  (setq g 0)
  (foreach w L
    (setq x (car w))
    (setq y (MATH::INT:func x))
    (setq g (+ g (* (cdr w) y)))
  )
  g
)
 
;;;=============================================================
;;; 获取埃米尔特系数                                            
;;;=============================================================
(defun Math::INT:GetHermite (n / xi wi EPS FI I ITS J M MAXIT P1 P2 P3 PIM4 PP Z Z1)
  (setq eps 1e-15)
  (setq PIM4 0.7511255444649425)
  (setq maxIt 10)
  (setq m (/ (1+ n) 2))
  (setq i 0)
  (while (< i m)
    (if (= i 0)
      (setq z (- (sqrt (+ n n 1)) (* 1.85575 (expt (+ n n 1) (/ -1 6.)))))
      (if (= i 1)
	(setq z (- z (/ (* 1.14 (expt n 0.426)) z)))
	(if (= i 2)
	  (setq z (- (* 1.86 z) (* 0.86 (cdr (assoc 0 xi)))))
	  (if (= i 3)
	    (setq z (- (* 1.91 z) (* 0.91 (cdr (assoc 1 xi)))))
	    (setq z (- (+ z z) (cdr (assoc (- i 2) xi))))
	  )
	)
      )
    )
    (setq its 0)
    (while (and (not (equal z z1 eps)) (< its MAXIT))
      (setq p1 pIM4)
      (setq p2 0.0)
      (setq j 0)
      (repeat n
	(setq p3 p2)
	(setq p2 p1)
	(setq p1 (- (* z p2 (sqrt (/ 2.0 (1+ j)))) (* p3 (sqrt (/ j (+ 1.0 j))))))
	(setq j (1+ j))
      )
      (setq pp (* p2 (sqrt (+ n n))))
      (setq z1 z)
      (setq z (- z1 (/ p1 pp)))
      (setq its (1+ its))
    )
    (setq fi (/ 2.0  pp pp))
    (setq xi (cons (cons i z) xi))
    (setq wi (cons (cons i fi) wi))
    (if (/= i (- n 1 i))
      (setq xi (cons (cons (- n 1 i) (- z)) xi)
	    wi (cons (cons (- n 1 i) fi) wi)
      )
    )
    (setq i (1+ i))
  )
  (MATH::INT:Bind xi wi)
)

方法三:高斯-埃米尔特积分

高斯-雅克比积分

;;;=============================================================
;;; 高斯-雅克比积分                                             
;;; 功能: 计算 f(x)*((1-x)^a)*((1+x)^b)的积分(在区间-1..1上)    
;;;=============================================================
(defun MATH::INT:Gauss-Jacobi (a b eps / ALF ARGS BET G N X Y flag g0)
  (if (setq args (UTI:InputBox))
    (progn
      (setq flag T)						;是否进行迭代
      (setq n 10)
      (while (and (< n 100) flag)	
        (setq alf (car args))
        (setq bet (cadr args))
        (setq g 0)
        (foreach w (MATH::INT:GetJacobiPolynomial n alf bet)
          (setq x (car w))
          (setq y (MATH::INT:func x))
          (setq g (+ g (* (cdr w) y)))
        )
	(if (equal g g0 eps)
	  (setq flag nil g0 g)
	  (setq n (+ n n) g0 g)
	)
      )
    )
  )
)
 
;;;=============================================================
;;; 高斯-雅克比积分                                             
;;; 功能: 计算高斯-雅克比积分系数项                             
;;;=============================================================
(defun MATH::INT:GetJacobiPolynomial (n alf bet /
				      xi wi A ALFBET AN B BN C EPS FI FLG I ITS
				      J MAXIT P1 P2 P3 PP R1 R2 R3 TEMP Z Z1)
  (setq maxIT 40)
  (setq eps 1e-15)
  (setq alf (float alf))
  (setq bet (float bet))
  (setq i 0)
  (repeat n									;for (i=0;i<n;i++) {
    (cond
      ( (= i 0)									;if (i == 0) {
        (setq an (/ alf n))							;an=alf/n;
        (setq bn (/ bet n))							;bn=bet/n;
        (setq r1 (* (1+ alf) (+ (/ 2.78 (+ 4.0 (* n n))) (/ (* 0.768 an) n))))	;r1=(1.0+alf)*(2.78/(4.0+n*n)+0.768*an/n);
        (setq r2 (+ 1.0 (* 1.48 n) (* 0.96 bn) (* 0.452 an an) (* 0.83 an bn)))	;r2=1.0+1.48*an+0.96*bn+0.452*an*an+0.83*an*bn;
        (setq z  (- 1 (/ r1 r2)))						;z=1.0-r1/r2;
      )
      ( (= i 1)									;else if (i == 1)
        (setq r1 (/ (+ 4.1 alf) (* (1+ alf) (1+ (* 0.156 alf)))))		;r1=(4.1+alf)/((1.0+alf)*(1.0+0.156*alf));
        (setq r2 (1+ (/ (* 0.06 (- n 8.0) (1+ (* 0.12 alf))) n)))		;r2=1.0+0.06*(n-8.0)*(1.0+0.12*alf)/n;
        (setq r3 (1+ (/ (* 0.012 bet (1+ (* 0.25 (abs alf)))) n)))		;r3=1.0+0.012*bet*(1.0+0.25*fabs(alf))/n;
        (setq z  (- z (* (- 1 z) r1 r2 r3)))					;z -= (1.0-z)*r1*r2*r3;
      )
      ( (= i 2)									;else if (i == 2)
        (setq r1 (/ (+ 1.67 (* 0.28 alf)) (1+ (* 0.37 alf))))			;r1=(1.67+0.28*alf)/(1.0+0.37*alf);
        (setq r2 (1+ (/ (* 0.22 (- n 8.0)) n)))					;r2=1.0+0.22*(n-8.0)/n;
        (setq r3 (1+ (/ (* 8.0 bet) (* (+ 6.28 bet) n n))))			;r3=1.0+8.0*bet/((6.28+bet)*n*n);
        (setq z  (- z (* (- (cdr (assoc 0 xi)) z) r1 r2 r3)))			;z -= (x[0]-z)*r1*r2*r3;
      )
      ( (= i (- n 2))								;else if (i == n-2)
        (setq r1 (/ (1+ (* 0.235 bet)) (+ 0.766 (* 0.119 bet))))		;r1=(1.0+0.37*bet)/(1.67+0.28*bet);
        (setq r2 (/ 1.0 (1+ (/ (* 0.639 (- n 4.0)) (1+ (* 0.71 (- n 4.0)))))))	;r2=1.0/(1.0+0.639*(n-4.0)/(1.0+0.71*(n-4.0)));
        (setq r3 (/ 1.0 (1+ (/ (* 20.0 alf) (* (+ 7.5 alf) n n))))) 		;r3=1.0/(1.0+20.0*alf/((7.5+alf)*n*n));
        (setq z  (+ z (* (- z (cdr (assoc (- n 4) xi))) r1 r2 r3))) 		;z += (z-x[n-4])*r1*r2*r3;
      )
      ( (= i (1- n))								;else if (i == n-1)
        (setq r1 (/ (1+ (* 0.37 bet))(+ 1.67 (* 0.28 bet))))			;r1=(1.0+0.37*bet)/(1.67+0.28*bet);
        (setq r2 (/ 1.0 (1+ (/ (* 0.22 (- n 8.0)) n))))				;r2=1.0/(1.0+0.22*(n-8.0)/n);
        (setq r3 (/ 1.0 (1+ (/ (* 8.0 alf) (* (+ 6.28 alf) n n)))))             ;r3=1.0/(1.0+8.0*alf/((6.28+alf)*n*n));
       	(setq z  (+ z (* (- z (cdr (assoc (- n 3) xi))) r1 r2 r3))) 		;z += (z-x[n-3])*r1*r2*r3;
      )
      (t									;else {
        (setq z (* 3 (- (cdr (assoc (1- i) xi)) (cdr (assoc (- i 2) xi)))))	;z=3.0*x[i-1]-3.0*x[i-2]+x[i-3];
        (setq z (+ z (cdr (assoc (- i 3) xi))))
      )
    )
    (setq alfbet (+ alf bet))							;alfbet=alf+bet;
    (setq its 1)
    (setq flg T)
    (while (and flg (<= its maxIt))						;for (its=1;its<=MAXIT;its++)
      (setq temp (+ 2.0 alfbet))						;temp=2.0+alfbet;
      (setq p1 (/ (+ (- alf bet) (* temp z)) 2.0))				;p1=(alf-bet+temp*z)/2.0;
      (setq p2 1.0)								;p2=1.0;
      (setq j 2)								
      (while (<= j n)								;for (j=2;j<=n;j++) {
	(setq p3 p2)								;p3=p2;
	(setq p2 p1)								;p2=p1;
	(setq temp (+ j j alfbet))						;temp=2*j+alfbet;
	(setq a (* 2 j (+ j alfbet) (- temp 2.0)))				;a=2*j*(j+alfbet)*(temp-2.0);
	(setq b (* (1- temp) (- (* alf alf) (* bet bet) (* temp (- 2 temp) z))));b=(temp-1.0)*(alf*alf-bet*bet+temp*(temp-2.0)*z);
	(setq c (* 2 (+ j -1 alf) (+ j -1 bet) temp))				;c=2.0*(j-1+alf)*(j-1+bet)*temp;
        (setq p1 (/ (- (* b p2) (* c p3)) a))					;p1=(b*p2-c*p3)/a;
	(setq j (1+ j))
      )
 
      (setq pp (/ (+ (* N (- ALF BET (* TEMP Z)) P1)(* 2 (+ N ALF)(+ N BET) P2))
		  (* TEMP (- 1.0 (* Z Z)))
	       )
      )										;pp=(n*(alf-bet-temp*z)*p1+2.0*(n+alf)*(n+bet)*p2)/(temp*(1.0-z*z));
      (setq z1 z)								;z1=z;
      (setq z (- z1 (/ p1 pp)))							;z=z1-p1/pp;
      (setq its (1+ its))
      (if (equal z z1 eps)
	(setq flg nil)
      )
    )
    (if (> its MAXIT)
      (princ "\nToo many iterations in gaujac!")				;if (its > MAXIT) nrerror("too many iterations in gaujac");
    )
    (setq xi (cons (cons i z) xi))						;x[i]=z;
    (setq fi (exp (- (Math::gammln (+ alf n))
		     (- (Math::gammln (+ bet n)))
		     (Math::gammln (1+ n))
		     (Math::gammln (+ n alfbet 1))
		  )
	     )
    )
    (setq fi (/ (* fi temp (expt 2.0 alfbet)) pp p2))				;w[i]=exp(gammln(alf+n)+gammln(bet+n)-gammln(n+1.0)-gammln(n+alfbet+1.0))*temp*pow(2.0,alfbet)/(pp*p2);	
    (setq wi (cons (cons i fi) wi))		
    (setq i (1+ i))
  )
  (MATH::INT:Bind xi wi)
)

四、高斯-拉盖尔积分
功能: 计算 e^(-x)*f(x)的广义积分(在区间0..INF上)

;;;=============================================================
;;; 高斯-拉盖尔积分                                             
;;; 功能: 计算 e^(-x)*f(x)的广义积分(在区间0..INF上)            
;;;=============================================================
(defun Math::INT:Gauss-Laguerre (a b eps / n flag L g h x y maxN)
  (setq n 2)
  (setq flag T)							;是否进行迭代
  (setq maxN 40)
  (while (and flag (< n maxN))	
    (setq g 0)
    (setq L (Math::INT:GetLaguerre n 0))
    (foreach w L
      (setq x (car w))
      (setq y (MATH::INT:func x))
      (setq g (+ g (* (cdr w) y)))
    )
    (if (equal g h eps)
      (setq flag nil)
      (setq n (+ n n)
	    h g
      )
    )
  )
  g
)
 
;;;=============================================================
;;; 获取拉盖尔系数                                              
;;;=============================================================
(defun Math::INT:GetLaguerre (n alf / xi wi AI EPS I ITS J MAXIT P1 P2 P3 PP X XI2 YI Z Z1)
  (setq maxIt 10)
  (setq eps 1e-15)
  (setq i 0)
  (repeat n
    (if	(= i 0)
      (setq z (/ (* (1+ alf) (+ 3 (* 0.92 alf)))
		 (+ 1 (* 2.4 n) (* 1.8 alf))
	      )
      )
      (if (= i 1)
	(setq
	  z (+ z (/ (+ 15 (* 6.25 alf)) (+ 1 (* 0.9 alf) (* 2.5 n))))
	)
	(setq ai  (1- i)
	      xi2 (cdr (assoc (- i 2) xi))
	      x	  (/ (*
		       (+
			 (/ (1+ (* 2.55 AI)) (* 1.9 AI))
			 (/ (* 1.26 AI ALF) (1+ (* 3.5 AI)))
		       )
		       (- Z XI2)
		     )
		     (1+ (* 0.3 ALF))
		  )
	      z	  (+ z x)
	)
      )
    )
    (setq its 0)
    (while (and (not (equal z z1 eps)) (< its maxIt))
      (setq p1 1)
      (setq p2 0)
      (setq j 0)
      (repeat n
	(setq p3 p2)
	(setq p2 p1)
	(setq p1 (/ (- (* (+ J J 1.0 ALF (- Z)) P2) (* (+ J ALF) P3))
		    (1+ J)
		 )
	)
	(setq j (1+ j))
      )
      (setq pp (/ (- (* n p1) (* p2 (+ n alf))) z))
      (setq z1 z)
      (setq z (- z1 (/ p1 pp)))
      (setq its (1+ its))
    )
    (setq xi (cons (cons i z) xi))
    (setq yi (-	(/ (exp (- (Math::gammln (+ alf n)) (Math::gammln n)))
		   (* pp n p2)
		)
	     )
    )
    (setq wi (cons (cons i yi) wi))
    (setq i (1+ i))
  )
  (MATH::INT:Bind xi wi)
)

五、高斯-切比雪夫积分
此方法用于针对 sqrt(1-x^2)*f(x)型的积分有很高的效率
;;;=============================================================
;;; 高斯-切比雪夫积分                                           
;;; 此方法用于针对 sqrt(1-x^2)*f(x)型的积分有很高的效率。       
;;; 输入: foo 函数名,arg 除自变量x外的参数列表, n迭代次数。     
;;; 输出: 此类积分的数值.                                       
;;; 说明: 此积分法效率较高,n取值一般10次左右就达到有效浮点精度.
;;; http://mathworld.wolfram.com/Chebyshev-GaussQuadrature.html 
;;;=============================================================
(defun Math::INT:Gauss-Chebyshev (a b eps / FI FLAG I N S0 SX WI XI)
  (setq n 3)							;一般来说迭代6-5次左右就可以达到浮点计算精度
  (setq flag T)							;是否进行迭代
  (while (and (< n 1000) flag)					;迭代次数不超过100
    (setq wi (/ pi n))						;所有的权值均为pi/n
    (setq sx 0)
    (setq i 1)
    (repeat n
      (setq xi (cos (/ (* pi (- i 0.5)) n)))			;x 项
      (setq fi (MATH::INT:func xi))				;x 项的函数值
      (setq sx (+ sx fi))
      (setq i  (1+ i))
    )
    (setq sx (* wi sx))						;统一乘以权值
    (if (equal sx s0 eps)					;是否满足精度要求
      (setq flag nil)						;是则不再迭代
      (setq n (+ n n)						;否则迭代次数倍增
	    s0 sx						;存储积分值
      )
    )
  )
  sx
)


下面是一般函数的积分方法:
六 、龙贝格积分法:
;;;=============================================================
;;; 龙贝格积分法                                                
;;; Romberg Integration                                         
;;; 输入: 函数名--foo (用符号表示,一般形式是 (foo x a b c ...) 
;;;       参数表--args ,除去自变量x的其它的参数列表            
;;;       下区间--Ra                                            
;;;       上区间--Rb                                            
;;;       容许误差--eps                                         
;;; 输出: 所求函数在区间段的积分                                
;;;=============================================================
(defun MATH::INT:Romberg (a b eps / EP H I K M N P Q S X Y Y0)
  (setq h (- b a))
  (setq y nil)
  (setq i 0)
  (repeat 20
    (setq y (cons (cons i 0.0) y))
    (setq i (1+ i))
  )
  (setq y (reverse y))
  (setq y0 (* h (+ (MATH::INT:func a) (MATH::INT:func b)) 0.5))
  (setq y (cons (cons 0 y0) (cdr y)))
  (setq	m  1
	n  1
	ep (1+ eps)
  )
  (while (and (>= ep eps) (<= m 19))
    (setq p 0.0)
    (setq i 0)
    (repeat n
      (setq x (+ a (* (+ i 0.5) h)))
      (setq p (+ p (MATH::INT:func x)))
      (setq i (1+ i))
    )
    (setq p (/ (+ (cdar y) (* h p)) 2.0))
    (setq s 1.0)
    (setq k 1)
    (repeat m
      (setq s (+ s s s s))
      (setq q (/ (- (* s p) (cdr (assoc (1- k) y))) (1- s)))
      (setq y (subst (cons (1- k) p) (assoc (1- k) y) y))
      (setq p q)
      (setq k (1+ k))
    )
    (setq ep (abs (- q (cdr (assoc (1- m) y)))))
    (setq m (1+ m))
    (setq y (subst (cons (1- m) q) (assoc (1- m) y) y))
    (setq n (+ n n))
    (setq h (/ h 2.0))
  )
  q
)

七 辛普森积分法。
;;;=============================================================
;;; 辛普森积分法                                                
;;; Simpson Integration                                         
;;;=============================================================
(defun MATH::INT:Simpson (a b eps / EP H ITER K N P S1 S2 T1 T2 X)
  (setq n 1)
  (setq h (- b a))
  (setq t1 (* h (+ (MATH::INT:func a) (MATH::INT:func b)) 0.5))
  (setq s1 t1)
  (setq ep (1+ eps))
  (setq iter 0)
  (while (and (>= ep eps) (< iter 50))
    (setq p 0.0)
    (setq k 0)
    (repeat n
      (setq x (+ a (* (+ k 0.5) h)))
      (setq p (+ p (MATH::INT:func x)))
      (setq k (1+ k))
    )
    (setq t2 (/ (+ t1 (* h p)) 2.))
    (setq s2 (/ (- (* 4.0 t2) t1) 3.))
    (setq ep (abs (- s2 s1)))
    (setq t1 t2)
    (setq s1 s2)
    (setq n (+ n n))
    (setq h (/ h 2))
    (setq iter (1+ iter))
  )
  s2
)

八、变步长梯形求积分法 
;;;=============================================================
;;; 变步长梯形求积分法                                          
;;; Trapezoidal Integration 1                                   
;;;=============================================================
(defun MATH::INT:Trapezia (a b eps / H K N P S T1 T2 X iter)
  (setq n 1)
  (setq h (- b a))
  (setq t1 (* h (+ (MATH::INT:func a) (MATH::INT:func b)) 0.5))
  (setq p (1+ eps))
  (setq iter 0)
  (while (and (>= p eps) (< iter 100))
    (setq s 0)
    (setq k 0)
    (repeat n
      (setq x (+ a (* (+ k 0.5) h)))
      (setq s (+ s (MATH::INT:func x)))
      (setq k (1+ k))
    )
    (setq t2 (/ (+ t1 (* h s)) 2.))
    (setq p (abs (- t1 t2)))
    (setq t1 t2)
    (setq n (+ n n))
    (setq h (/ h 2))
    (setq iter (1+ iter))
  )
  t2
)


九、步长积分法。
;;;=============================================================
;;; 步长积分法                                                  
;;; Trapezoidal Integration 2                                   
;;;=============================================================
(defun MATH::INT:Trapzd (a b n / DEL IT SUM TNM X s)
  (if (= n 1)
    (setq s (* 0.5 (- b a) (+ (MATH::INT:func a) (MATH::INT:func b))))
    (progn
      (setq it 1)
      (repeat (- n 2)
	(setq it (lsh it 1))
      )
      (setq tnm it)
      (setq del (/ (- b a) tnm))
      (setq x (+ a (* 0.5 del)))
      (setq sum 0.0)
      (repeat it
	(setq sum (+ sum (MATH::INT:func x)))
	(setq x (+ x del))
      )
      (setq s (* 0.5 (+ s (/ (* (- b a) sum) tnm))))
    )
  )
)


十、自适应积分法
;;;=============================================================
;;; 自适应求积分法                                              
;;; Self-adapting Trapezia Integration                          
;;;=============================================================
(defun MATH::INT:Atrapezia (a b eps / F0 F1 H S T0 TT d)
  (setq d 1e-4)
  (setq h (- b a))
  (setq TT '(0. . 0.))
  (setq f0 (MATH::INT:func a))
  (setq f1 (MATH::INT:func b))
  (setq t0 (* h (+ f0 f1) 0.5))
  (car (MATH::INT:ppp a b h f0 f1 t0 eps d tt))
)
 
(defun MATH::INT:PPP (x0 x1 h f0 f1 t0 eps d tt / EPS1 F G P T1 T2 T3 X X2)
  (setq x (+ x0 (* h 0.5)))
  (setq f (MATH::INT:func x))
  (setq t1 (* h (+ f0 f) 0.25))
  (setq t2 (* h (+ f1 f) 0.25))
  (setq p (abs (- t0 t1 t2)))
  (if (or (< p eps) (< (* 0.5 h) d))
    (cons (+ (car tt) t1 t2) (cdr tt))
    (progn
      (setq g (* h 0.5))
      (setq eps1 (/ eps 1.4))
      (setq t3 (MATH::INT:ppp x0 x g f0 f t1 eps1 d tt))
      (setq t3 (MATH::INT:ppp x x1 g f f1 t2 eps1 d t3))
    )
  )
)

一些相关函数:
;;;=============================================================
;;; 功能: 第一类椭圆积分                                        
;;; 输入: Phi 和 k < 1                                          
;;; 输出: 所求椭圆积分值                                        
;;;=============================================================
(defun Math:Elliptic_Integral_1 (phi kCoff / )
  (defun MATH::INT:func (x / s)
    (setq s (* kCoff (sin x)))
    (/ 1.0 (sqrt (* (- 1 s) (1+ s))))
  )
  (MATH::INT:Romberg 0 phi 1e-15)
)
 
;;;=============================================================
;;; 功能: 第二类椭圆积分                                        
;;; 输入: Phi 和 k < 1                                          
;;; 输出: 所求椭圆积分值                                        
;;;=============================================================
(defun Math:Elliptic_Integral_2 (phi kCoff / )
  (defun MATH::INT:func (x /)
    (setq x (* kCoff (sin x)))
    (sqrt (* (- 1 x) (1+ x)))
  )
  (MATH::INT:Romberg 0 phi 1e-15)
)
 
;;;=============================================================
;;; 阶乘                                                        
;;;=============================================================
(defun Math::frac (n)
  (if (= n 0)
    1.
    (if	(= n 1)
      1.
      (* n (Math::frac (1- n)))
    )
  )
)
 
;;;=============================================================
;;; 伽玛函数                                                    
;;;=============================================================
(defun Math::gammln (xx / x y tmp ser j cof)
  (setq	cof '(76.18009172947146		 -86.50532032941677
	      24.01409824083091		 -1.231739572450155
	      0.1208650973866179e-2	 -0.5395239384953e-5
	     )
  )
  (setq x xx)
  (setq y x)
  (setq tmp (+ x 5.5))
  (setq tmp (- tmp (* (+ 0.5 x) (log tmp))))
  (setq ser 1.000000000190015)
  (setq j 0)
  (repeat 6
    (setq y (1+ y))
    (setq ser (+ ser (/ (nth j cof) y)))
    (setq j (1+ j))
  )
  (- (log (/ (* 2.5066282746310005 ser) x)) tmp)
)
 
;;;=============================================================
;;; 次函数用于比较排序                                          
;;;=============================================================
(defun MATH::INT:funcSort (e1 e2)
  (< (car e1) (car e2))
)
 
;;;=============================================================
;;; 组合两个系数                                                
;;;=============================================================
(defun MATH::INT:Bind (xi wi)
  (setq xi (vl-sort xi 'MATH::INT:funcSort))
  (setq wi (vl-sort wi 'MATH::INT:funcSort))
  (setq xi (mapcar 'cdr xi))
  (setq wi (mapcar 'cdr wi))
  (mapcar 'cons xi wi)
)



下面是用对话框创建的求积分法的例程:

;;;=============================================================
;;; 用各种方法求积分的程序                                      
;;;=============================================================
(defun C:Quadrature (/ ID OK DCL_FILE)
  (setq id (load_dialog (setq Dcl_File (MATH::INT:Write_Dcl))))	;从对话框中得到表达式
  (vl-file-delete Dcl_File)					;删除临时对话框文件
  (setq ok 2)
  (if (new_dialog "dcl_Integration" id)
    (progn
      (VL-CATCH-ALL-APPLY 'MATH::INT:GetSettings)		;读取默认数据
      (action_tile "help" "(MATH::INT:Help 1)")			;帮助
      (foreach k '(0 1 2 3 4 5 6 7 8)
	(setq k (strcat "K" (itoa k)))
        (action_tile k "(MATH::INT:OnBtn $key)")                ;按钮动作,对应相对的积分方法
      )		   
      (setq ok (start_dialog))
    )
  )
  (unload_dialog ID)
  (princ)
)
 
(defun C:JF (/)
  (VL-CATCH-ALL-APPLY 'C:Quadrature)
  (princ)
)
 
;;;=============================================================
;;; 从环境变量读取上次数据                                      
;;;=============================================================
(defun MATH::INT:GetSettings (/ data)
  (if (setq Data (getenv "Intergration"))
    (foreach k (read data)
      (set_tile (car k) (cdr k))
    )
  )
)
 
;;;=============================================================
;;; 检查对话框输入                                              
;;;=============================================================
(defun MATH::INT:CheckInput (symS symA symB symN / e f)
  (setq e (exp 1))
  (set symS (get_tile "F"))
  (set symA (MATH::INT:MyRead (get_tile "A")))
  (set symB (MATH::INT:MyRead (get_tile "B")))
  (set symN (MATH::INT:MyRead (get_tile "N")))
  (setq f (CAL:Expr2Func (eval s) 'MATH::INT:func '(x)))
  (apply 'and (mapcar 'eval '(F symS symA symB symN)))   
)
 
 
;;;=============================================================
;;; 按钮动作,对应相应的函数求积分                              
;;;=============================================================
(defun MATH::INT:OnBtn (key / DATA s N a b OldZIN m EPS RET tm0 map msg foo tmp idx)
  (setq m (VL-CATCH-ALL-APPLY 'MATH::INT:CheckInput '(s a b n)))
  (if (or (vl-catch-all-error-p m) (not m) (equal a b 1e-8))
    (if (vl-catch-all-error-p m) 
      (set_tile "info" (vl-catch-all-error-message m))
      (set_tile "info" "无效输入!")
    )
    (progn
      ;;如果精度过高,设置为15位的精度
      (if (> n 20)
	(setq n 15)
	(setq n (fix (abs n)))
      )
      ;;如果上区间小于下区间,则交换区间
      (if (< b a)
	(setq tmp a a b b tmp)
      )
      ;;记住对话框输入,用于下次
      (setq OldZIN (getvar "DIMZIN"))
      (setvar "DIMZIN" 8)
      (set_tile "N" (itoa n))
      (set_tile "A" (rtos a 2 20))
      (set_tile "B" (rtos b 2 20))
      (setq data (list (cons "F" s)
		       (cons "N" (itoa n))
		       (cons "A" (rtos a 2 20))
		       (cons "B" (rtos b 2 20))
		 )
      )
      (setvar "DIMZIN" OldZIN)
      (setenv "Intergration" (VL-PRIN1-TO-STRING data))
      ;;开始计算积分
      (setq eps (expt 0.1 n))
      (setq tm0 (getvar "TDUSRTIMER"))
      (setq map (MATH::INT:GetMethods))				;积分计算方法集
      (setq idx (atoi (substr key 2)))				
      (setq foo (nth idx map))					;获取计算积分的函数
      (setq ret (VL-CATCH-ALL-APPLY foo (list a b eps)))        ;获取积分值
      (if (vl-catch-all-error-p ret)
	(set_tile "info" (vl-catch-all-error-message ret))	;求解过程发生了错误
        (if (null ret)
	  (set_tile "info" "发生了错误,求值结果为空!")	
	  (progn
	    ;(MATH::INT:Bench 100 a b eps)
            (setq ret (rtos ret 2 20))
            (setq msg (get_attr key "label"))
            (setq msg (strcat msg "求的结果为:" ret))
            (set_tile (strcat "R" (itoa idx)) ret)             	;显示求解结果
            (set_tile "info" msg)				;显示求解结果
            (princ (strcat "\n" msg))				;打印求解结果
            (princ "\n费时:")
            (princ (* (- (getvar "TDUSRTIMER") tm0) 86400))
            (princ "秒.")
	  )
        )
      )
    )
  )
)
 
;;;=============================================================
;;; 各种积分测速                                                
;;;=============================================================
(defun MATH::INT:Bench (n a b eps)
  (UTI:BENCH
    n
    (list
      (list 'MATH::INT:Romberg a b eps)
      (list 'MATH::INT:Gauss-Legendre a b eps)
      (list 'MATH::INT:Simpson a b eps)
    )
  )
)
 
;;;=============================================================
;;; 积分方法集                                                  
;;;=============================================================
(defun MATH::INT:GetMethods ()
  '(MATH::INT:Romberg
    MATH::INT:Simpson
    MATH::INT:Atrapezia
    MATH::INT:Trapezia
    MATH::INT:Gauss-Legendre
    Math::INT:Gauss-Chebyshev
    Math::INT:Gauss-Laguerre
    Math::INT:Gauss-Hermite
    MATH::INT:Gauss-Jacobi
   )
)
 
;;;=============================================================
;;; 表达式求值,也可以用cal函数                                 
;;;=============================================================
(defun MATH::INT:MyRead (str / e)
  (setq e (exp 1))
  (CAL:Expr2Value str)
)
 
;;;=============================================================
;;; 帮助和说明: help and instruction                            
;;;=============================================================
(defun MATH::INT:Help (n)
  (if (= n 1)
    (if	(= "CHS" (getvar "Locale"))
      (alert
	"函数式只接受符号x为变量,不规范很可能出错!
	\n函数可以LISP内置的数学函数,也可以自定义函数!
	\n指数用^表示,+-*/表示加减乘除,乘号不能省略。
	\n程序能采用多种方法求积,一般来说龙贝格积分法最快。
	\n有什么问题email: highflybird@qq.com
	\n作者: highflybird 日期2019.07"
      )
      (alert
	"Standard expression only accepts \"x\" as a variale!
	\nThe fastest is Romberg Integration,the slowest is Trapezoidal rule(Be careful!).
	\nRecommendation:Don't set a high precision at first,promote it step by step.
	\nEspically for the Trapezoidal rule, It won't work well on some circumstances.
	\nIt's an Open Source Software. Thanks for your advice or bug reports.
	\nAuthor: highflybird  Email: highflybird@qq.com  Date:2019.07."
      )
    )
    (set_tile "info" "表达式非法或者空输入.")
  )
)

对话框的制作:
;;;=============================================================
;;; 输入对话框                                                  
;;;=============================================================
(defun UTI:Inputbox (/ str wcs ret)
  (setq	str "Function GetNumbers()
  	     GetNumbers=inputbox(\"请输入两个参数,中间用空格隔开:\",\"输入框\")
             End Function"
  )
  (if
    (or
      (setq wcs (vlax-create-object "Aec32BitAppServer.AecScriptControl.1"))
      (setq wcs (vlax-create-object "ScriptControl"))
    )
    (progn 
      (vlax-put-property wcs "language" "VBScript")
      (vlax-invoke wcs 'addcode str)
      (if (setq ret (vlax-invoke wcs 'run "GetNumbers"))
	(setq ret (strcat "(" ret ")")
	      ret (read ret)
	)
      )
      (vlax-release-object wcs)
      (if
	(and
	  (= 2 (length ret))
	  (or (= 'INT (type (car ret))) (= 'REAL (type (car ret))))
	  (or (= 'INT (type (cadr ret))) (= 'REAL (type (cadr ret))))
	)
	ret
      )
    )
  )
)
 
;;;=============================================================
;;; 写对话框到文件用于程序                                      
;;;=============================================================
(defun MATH::INT:Write_Dcl (/ Dcl_File file str)
  (setq Dcl_File (vl-filename-mktemp nil nil ".Dcl"))
  (setq file (open Dcl_File "w"))
  (princ
    "dcl_Integration : dialog {
	label = \"数值积分LISP版  v1.2\";
	: boxed_column {
          width = 60;
	  fixed_width = true;
	  : edit_box {
	    key=\"F\";
	    label= \"函数:\";
	  }
	  : row {
	    : edit_box {
	      key=\"A\";
	      label= \"下限:\";
	    }
	    : edit_box {
	      key=\"B\";
	      label= \"上限:\";
	    }      
	    : edit_box {
	      key=\"N\";
	      label = \"精确位数:\";
	      value = 8;
	      edit_width = 2;
	      fixed_width = true;
	    }
	  }
	  spacer_1;
	}
	: row {
	  : boxed_column {
	    label = \"计算方法:\";
	    : button {
	      key = \"K0\";
	      label = \"龙贝格积分法\";
	    }
	    : button {
	      key = \"K1\";
	      label = \"辛普森积分法\";
	    }
	    : button {
	      key = \"K2\";
	      label = \"自适应积分法\";
	    }
	    : button {
	      key = \"K3\";
	      label = \"变步长梯形法\";
	    }
	    : button {
	      key = \"K4\";
	      label = \"高斯-勒让德法\";
	    }
	    : button {
	      key = \"K5\";
	      label = \"高斯-切比雪夫法\";
	    }
	    : button {
	      key = \"K6\";
	      label = \"高斯-拉盖尔法\";
	    }
	    : button {
	      key = \"K7\";
	      label = \"高斯-埃尔米特法\";
	    }
	    : button {
	      key = \"K8\";
	      label = \"高斯-雅克比法\";
	    }
	    spacer;
	  }
	  : boxed_column {
	    width = 32;
	    fixed_width = true;
	    label = \"计算结果:\";
	    : text {
	      key = \"R0\";
	    }
	    : text {
	      key = \"R1\";
	    }
	    : text {
	      key = \"R2\";
	    }
	    : text {
	      key = \"R3\";
	    }
	    : text {
	      key = \"R4\";
	    }
	    : text {
	      key = \"R5\";
	    }
	    : text {
	      key = \"R6\";
	    }
	    : text {
	      key = \"R7\";
	    }
	    : text {
	      key = \"R8\";
	    }
	    spacer;
	  }
	}
	ok_cancel_help;
	//ok_cancel_help_errtile;
	: text {
	  key = \"info\";
	  label = \"Copyright \\u+00A9 2007-2019 Highflybird. All rights reserved.\";
	  width = 20;
	}
    }  "
    file
  )
  (close file)
  Dcl_File
)
 
;;;=============================================================
;;; 下面的就不用介绍了                                          
;;;=============================================================
(vl-load-com)
(if (= "CHS" (getvar "Locale"))
  (prompt "输入命令: JF")
  (prompt "Please enter: Quadrature")
)
(c:JF)
(princ)


一些应用,如下面网友的提问并解答:

http://bbs.mjtd.com/forum.php?mod=viewthread&tid=179809&extra=&highlight=%BB%FD%B7%D6&page=1

椭球体球缺的表面积计算
积分公式用LISP程序表达并计算出来
这是一个网上的
有结果和公式

用上面的介绍的一些函数可以得出其面积和体积:

用LISP编写了一个求积分的程序:

里面采用了各种方法求积分和各种类型的积分。下面我把各种方法的源码贴出。

方法一: 勒让德-高斯积分法。

;;;=============================================================
;;; 计算勒让德-高斯求积函数的系数 ,如下面的6次项系数           
;;; p0(x) = 1                                                   
;;; p1(x) = x                                                   
;;; p2(x) = (3*x^2-1)/2                                         
;;; p3(x) = (5*x^3-3*x)/2                                       
;;; p4(x) = (35*x^4-30*x^2+3)/8                                 
;;; p5(x) = (63*x^5-70*x^3+15*x)/8                              
;;; p6(x) = (231*x^6-315*x^4+105*x^2-5)/16                      
;;; x-2                                                         
;;; x-1                                                         
;;; 9x-5                                                        
;;; 216*x^2-216*x+49                                            
;;; 45000*x^2-32200*x+5103                                      
;;; 2025000*x^3-2025000*x^2+629325*x-58564                      
;;; 142943535000*x^3-1130712534000*x^2+27510743799*x-1976763932 
;;; 输入: x1,x2 区间(一般来说是-1..1),n 迭代次数,eps迭代精度   
;;; 输出: 勒让德-高斯求积函数的系数,用点表集表示                
;;; http://mathworld.wolfram.com/Legendre-GaussQuadrature.html  
;;;=============================================================
(defun Math::Int:Legendre_Polynomial (x1 x2 n eps / xi wi FI I ITER J M P1 P2 P3 PP XL XM Z Z1)
  (setq m (/ (1+ n) 2))
  (setq xm (* 0.5 (+ x2 x1)))
  (setq xl (* 0.5 (- x2 x1)))
  (setq i 1)
  (repeat m
    (setq z (cos (/ (* pi (- i 0.25)) (+ n 0.5))))
    (setq iter 0)
    (while (and (not (equal z z1 eps)) (< iter 1000))
      (setq p1 1.0)
      (setq p2 0.0)
      (setq j 1)
      (repeat n
	(setq p3 p2)
	(setq p2 p1)
	(setq p1 (/ (- (* z p2 (+ j j -1.)) (* (1- j) p3)) j))
	(setq j (1+ j))
      )
      (setq pp (/ (* n (- (* z p1) p2)) (1- (* z z))))
      (setq z1 z)
      (setq z (- z1 (/ p1 pp)))
      (setq iter (1+ iter))
    )
    (setq fi (/ xl 0.5 (- 1 (* z z)) pp pp))
    (setq xi (cons (cons (1- i) (- xm (* xl z))) xi))
    (setq wi (cons (cons (1- i) fi) wi))
    (if (/= (1- i) (- n i))
      (setq xi (cons (cons (- n i) (+ xm (* xl z))) xi)
	    wi (cons (cons (- n i) fi) wi)
      )
    )
    (setq i (1+ i)) 
  )
  (MATH::INT:Bind xi wi)
)
 
;;;=============================================================
;;; 勒让德-高斯求积函数                                         
;;;=============================================================
(defun MATH::INT:Gauss-Legendre (a b eps / AA BB EP FX G H I L M P S W X)
  (setq l '((-0.93246951420315202787 . 0.17132449237917034504)  
	    (-0.66120938646626451363 . 0.36076157304813860756)     
	    (-0.23861918608319690859 . 0.46791393457269104739)  
	    ( 0.23861918608319690859 . 0.46791393457269104739)     
            ( 0.66120938646626451363 . 0.36076157304813860756)     
	    ( 0.93246951420315202787 . 0.17132449237917034504)
	  )
  )                                                         
  ;(setq L (Math::Int:Legendre_Polynomial -1 1 100 2e-20))
  (setq m 1)
  (setq h (- b a))
  (setq s (abs (* 0.001 h)))
  (setq p 1e100) 	
  (setq ep (1+ eps))
  (while (and (>= ep eps) (> (abs h) s))
    (setq g 0)
    (setq i 1)
    (repeat m
      (setq bb (+ a (* i h)))
      (setq aa (- bb h))
      (setq w 0)
      (foreach k l
	(setq x (* 0.5 (+ bb aa (* (- bb aa) (car k)))))
	(setq fx (MATH::INT:func x))
	(setq w (+ w (* fx (cdr k))))
      )
      (setq g (+ g w))
      (setq i (1+ i))
    ) 
    (setq g (* g h 0.5))
    (setq ep (/ (abs (- g p)) (1+ (abs g))))
    (setq p g)
    (setq m (1+ m))
    (setq h (/ (- b a) m))
  )
  g
)
 
;;;=============================================================
;;; 勒让德-高斯求积函数(另一方法,慢些)                         
;;;=============================================================
(defun MATH::INT:Gauss-Legendre1 (a b eps / FLAG G H L N X Y)
  (setq n 1)                                                      
  (setq flag T)							;是否进行迭代
  (while (and (< n 100) flag)	
    (setq g 0)
    (setq L (Math::Int:Legendre_Polynomial a b n eps))
    (foreach w L
      (setq x (car w))
      (setq y (MATH::INT:func x))
      (setq g (+ g (* (cdr w) y)))
    )
    (if (equal g h eps)
      (setq flag nil)
      (setq n (+ n n)
	    h g
      )
    )
  )
  g
)

方法二:高斯-埃米尔特积分

;;;=============================================================
;;; 高斯-埃米尔特积分                                           
;;; 功能: 计算 e^(-x^2)*f(x)的广义积分(在区间-INF..INF上)       
;;;=============================================================
(defun Math::INT:Gauss-Hermite (a b eps / n L g x y)
  (setq n 100)
  (setq L (Math::INT:GetHermite n))
  (setq g 0)
  (foreach w L
    (setq x (car w))
    (setq y (MATH::INT:func x))
    (setq g (+ g (* (cdr w) y)))
  )
  g
)
 
;;;=============================================================
;;; 获取埃米尔特系数                                            
;;;=============================================================
(defun Math::INT:GetHermite (n / xi wi EPS FI I ITS J M MAXIT P1 P2 P3 PIM4 PP Z Z1)
  (setq eps 1e-15)
  (setq PIM4 0.7511255444649425)
  (setq maxIt 10)
  (setq m (/ (1+ n) 2))
  (setq i 0)
  (while (< i m)
    (if (= i 0)
      (setq z (- (sqrt (+ n n 1)) (* 1.85575 (expt (+ n n 1) (/ -1 6.)))))
      (if (= i 1)
	(setq z (- z (/ (* 1.14 (expt n 0.426)) z)))
	(if (= i 2)
	  (setq z (- (* 1.86 z) (* 0.86 (cdr (assoc 0 xi)))))
	  (if (= i 3)
	    (setq z (- (* 1.91 z) (* 0.91 (cdr (assoc 1 xi)))))
	    (setq z (- (+ z z) (cdr (assoc (- i 2) xi))))
	  )
	)
      )
    )
    (setq its 0)
    (while (and (not (equal z z1 eps)) (< its MAXIT))
      (setq p1 pIM4)
      (setq p2 0.0)
      (setq j 0)
      (repeat n
	(setq p3 p2)
	(setq p2 p1)
	(setq p1 (- (* z p2 (sqrt (/ 2.0 (1+ j)))) (* p3 (sqrt (/ j (+ 1.0 j))))))
	(setq j (1+ j))
      )
      (setq pp (* p2 (sqrt (+ n n))))
      (setq z1 z)
      (setq z (- z1 (/ p1 pp)))
      (setq its (1+ its))
    )
    (setq fi (/ 2.0  pp pp))
    (setq xi (cons (cons i z) xi))
    (setq wi (cons (cons i fi) wi))
    (if (/= i (- n 1 i))
      (setq xi (cons (cons (- n 1 i) (- z)) xi)
	    wi (cons (cons (- n 1 i) fi) wi)
      )
    )
    (setq i (1+ i))
  )
  (MATH::INT:Bind xi wi)
)

方法三:高斯-埃米尔特积分

高斯-雅克比积分

;;;=============================================================
;;; 高斯-雅克比积分                                             
;;; 功能: 计算 f(x)*((1-x)^a)*((1+x)^b)的积分(在区间-1..1上)    
;;;=============================================================
(defun MATH::INT:Gauss-Jacobi (a b eps / ALF ARGS BET G N X Y flag g0)
  (if (setq args (UTI:InputBox))
    (progn
      (setq flag T)						;是否进行迭代
      (setq n 10)
      (while (and (< n 100) flag)	
        (setq alf (car args))
        (setq bet (cadr args))
        (setq g 0)
        (foreach w (MATH::INT:GetJacobiPolynomial n alf bet)
          (setq x (car w))
          (setq y (MATH::INT:func x))
          (setq g (+ g (* (cdr w) y)))
        )
	(if (equal g g0 eps)
	  (setq flag nil g0 g)
	  (setq n (+ n n) g0 g)
	)
      )
    )
  )
)
 
;;;=============================================================
;;; 高斯-雅克比积分                                             
;;; 功能: 计算高斯-雅克比积分系数项                             
;;;=============================================================
(defun MATH::INT:GetJacobiPolynomial (n alf bet /
				      xi wi A ALFBET AN B BN C EPS FI FLG I ITS
				      J MAXIT P1 P2 P3 PP R1 R2 R3 TEMP Z Z1)
  (setq maxIT 40)
  (setq eps 1e-15)
  (setq alf (float alf))
  (setq bet (float bet))
  (setq i 0)
  (repeat n									;for (i=0;i<n;i++) {
    (cond
      ( (= i 0)									;if (i == 0) {
        (setq an (/ alf n))							;an=alf/n;
        (setq bn (/ bet n))							;bn=bet/n;
        (setq r1 (* (1+ alf) (+ (/ 2.78 (+ 4.0 (* n n))) (/ (* 0.768 an) n))))	;r1=(1.0+alf)*(2.78/(4.0+n*n)+0.768*an/n);
        (setq r2 (+ 1.0 (* 1.48 n) (* 0.96 bn) (* 0.452 an an) (* 0.83 an bn)))	;r2=1.0+1.48*an+0.96*bn+0.452*an*an+0.83*an*bn;
        (setq z  (- 1 (/ r1 r2)))						;z=1.0-r1/r2;
      )
      ( (= i 1)									;else if (i == 1)
        (setq r1 (/ (+ 4.1 alf) (* (1+ alf) (1+ (* 0.156 alf)))))		;r1=(4.1+alf)/((1.0+alf)*(1.0+0.156*alf));
        (setq r2 (1+ (/ (* 0.06 (- n 8.0) (1+ (* 0.12 alf))) n)))		;r2=1.0+0.06*(n-8.0)*(1.0+0.12*alf)/n;
        (setq r3 (1+ (/ (* 0.012 bet (1+ (* 0.25 (abs alf)))) n)))		;r3=1.0+0.012*bet*(1.0+0.25*fabs(alf))/n;
        (setq z  (- z (* (- 1 z) r1 r2 r3)))					;z -= (1.0-z)*r1*r2*r3;
      )
      ( (= i 2)									;else if (i == 2)
        (setq r1 (/ (+ 1.67 (* 0.28 alf)) (1+ (* 0.37 alf))))			;r1=(1.67+0.28*alf)/(1.0+0.37*alf);
        (setq r2 (1+ (/ (* 0.22 (- n 8.0)) n)))					;r2=1.0+0.22*(n-8.0)/n;
        (setq r3 (1+ (/ (* 8.0 bet) (* (+ 6.28 bet) n n))))			;r3=1.0+8.0*bet/((6.28+bet)*n*n);
        (setq z  (- z (* (- (cdr (assoc 0 xi)) z) r1 r2 r3)))			;z -= (x[0]-z)*r1*r2*r3;
      )
      ( (= i (- n 2))								;else if (i == n-2)
        (setq r1 (/ (1+ (* 0.235 bet)) (+ 0.766 (* 0.119 bet))))		;r1=(1.0+0.37*bet)/(1.67+0.28*bet);
        (setq r2 (/ 1.0 (1+ (/ (* 0.639 (- n 4.0)) (1+ (* 0.71 (- n 4.0)))))))	;r2=1.0/(1.0+0.639*(n-4.0)/(1.0+0.71*(n-4.0)));
        (setq r3 (/ 1.0 (1+ (/ (* 20.0 alf) (* (+ 7.5 alf) n n))))) 		;r3=1.0/(1.0+20.0*alf/((7.5+alf)*n*n));
        (setq z  (+ z (* (- z (cdr (assoc (- n 4) xi))) r1 r2 r3))) 		;z += (z-x[n-4])*r1*r2*r3;
      )
      ( (= i (1- n))								;else if (i == n-1)
        (setq r1 (/ (1+ (* 0.37 bet))(+ 1.67 (* 0.28 bet))))			;r1=(1.0+0.37*bet)/(1.67+0.28*bet);
        (setq r2 (/ 1.0 (1+ (/ (* 0.22 (- n 8.0)) n))))				;r2=1.0/(1.0+0.22*(n-8.0)/n);
        (setq r3 (/ 1.0 (1+ (/ (* 8.0 alf) (* (+ 6.28 alf) n n)))))             ;r3=1.0/(1.0+8.0*alf/((6.28+alf)*n*n));
       	(setq z  (+ z (* (- z (cdr (assoc (- n 3) xi))) r1 r2 r3))) 		;z += (z-x[n-3])*r1*r2*r3;
      )
      (t									;else {
        (setq z (* 3 (- (cdr (assoc (1- i) xi)) (cdr (assoc (- i 2) xi)))))	;z=3.0*x[i-1]-3.0*x[i-2]+x[i-3];
        (setq z (+ z (cdr (assoc (- i 3) xi))))
      )
    )
    (setq alfbet (+ alf bet))							;alfbet=alf+bet;
    (setq its 1)
    (setq flg T)
    (while (and flg (<= its maxIt))						;for (its=1;its<=MAXIT;its++)
      (setq temp (+ 2.0 alfbet))						;temp=2.0+alfbet;
      (setq p1 (/ (+ (- alf bet) (* temp z)) 2.0))				;p1=(alf-bet+temp*z)/2.0;
      (setq p2 1.0)								;p2=1.0;
      (setq j 2)								
      (while (<= j n)								;for (j=2;j<=n;j++) {
	(setq p3 p2)								;p3=p2;
	(setq p2 p1)								;p2=p1;
	(setq temp (+ j j alfbet))						;temp=2*j+alfbet;
	(setq a (* 2 j (+ j alfbet) (- temp 2.0)))				;a=2*j*(j+alfbet)*(temp-2.0);
	(setq b (* (1- temp) (- (* alf alf) (* bet bet) (* temp (- 2 temp) z))));b=(temp-1.0)*(alf*alf-bet*bet+temp*(temp-2.0)*z);
	(setq c (* 2 (+ j -1 alf) (+ j -1 bet) temp))				;c=2.0*(j-1+alf)*(j-1+bet)*temp;
        (setq p1 (/ (- (* b p2) (* c p3)) a))					;p1=(b*p2-c*p3)/a;
	(setq j (1+ j))
      )
 
      (setq pp (/ (+ (* N (- ALF BET (* TEMP Z)) P1)(* 2 (+ N ALF)(+ N BET) P2))
		  (* TEMP (- 1.0 (* Z Z)))
	       )
      )										;pp=(n*(alf-bet-temp*z)*p1+2.0*(n+alf)*(n+bet)*p2)/(temp*(1.0-z*z));
      (setq z1 z)								;z1=z;
      (setq z (- z1 (/ p1 pp)))							;z=z1-p1/pp;
      (setq its (1+ its))
      (if (equal z z1 eps)
	(setq flg nil)
      )
    )
    (if (> its MAXIT)
      (princ "\nToo many iterations in gaujac!")				;if (its > MAXIT) nrerror("too many iterations in gaujac");
    )
    (setq xi (cons (cons i z) xi))						;x[i]=z;
    (setq fi (exp (- (Math::gammln (+ alf n))
		     (- (Math::gammln (+ bet n)))
		     (Math::gammln (1+ n))
		     (Math::gammln (+ n alfbet 1))
		  )
	     )
    )
    (setq fi (/ (* fi temp (expt 2.0 alfbet)) pp p2))				;w[i]=exp(gammln(alf+n)+gammln(bet+n)-gammln(n+1.0)-gammln(n+alfbet+1.0))*temp*pow(2.0,alfbet)/(pp*p2);	
    (setq wi (cons (cons i fi) wi))		
    (setq i (1+ i))
  )
  (MATH::INT:Bind xi wi)
)

四、高斯-拉盖尔积分
功能: 计算 e^(-x)*f(x)的广义积分(在区间0..INF上)

;;;=============================================================
;;; 高斯-拉盖尔积分                                             
;;; 功能: 计算 e^(-x)*f(x)的广义积分(在区间0..INF上)            
;;;=============================================================
(defun Math::INT:Gauss-Laguerre (a b eps / n flag L g h x y maxN)
  (setq n 2)
  (setq flag T)							;是否进行迭代
  (setq maxN 40)
  (while (and flag (< n maxN))	
    (setq g 0)
    (setq L (Math::INT:GetLaguerre n 0))
    (foreach w L
      (setq x (car w))
      (setq y (MATH::INT:func x))
      (setq g (+ g (* (cdr w) y)))
    )
    (if (equal g h eps)
      (setq flag nil)
      (setq n (+ n n)
	    h g
      )
    )
  )
  g
)
 
;;;=============================================================
;;; 获取拉盖尔系数                                              
;;;=============================================================
(defun Math::INT:GetLaguerre (n alf / xi wi AI EPS I ITS J MAXIT P1 P2 P3 PP X XI2 YI Z Z1)
  (setq maxIt 10)
  (setq eps 1e-15)
  (setq i 0)
  (repeat n
    (if	(= i 0)
      (setq z (/ (* (1+ alf) (+ 3 (* 0.92 alf)))
		 (+ 1 (* 2.4 n) (* 1.8 alf))
	      )
      )
      (if (= i 1)
	(setq
	  z (+ z (/ (+ 15 (* 6.25 alf)) (+ 1 (* 0.9 alf) (* 2.5 n))))
	)
	(setq ai  (1- i)
	      xi2 (cdr (assoc (- i 2) xi))
	      x	  (/ (*
		       (+
			 (/ (1+ (* 2.55 AI)) (* 1.9 AI))
			 (/ (* 1.26 AI ALF) (1+ (* 3.5 AI)))
		       )
		       (- Z XI2)
		     )
		     (1+ (* 0.3 ALF))
		  )
	      z	  (+ z x)
	)
      )
    )
    (setq its 0)
    (while (and (not (equal z z1 eps)) (< its maxIt))
      (setq p1 1)
      (setq p2 0)
      (setq j 0)
      (repeat n
	(setq p3 p2)
	(setq p2 p1)
	(setq p1 (/ (- (* (+ J J 1.0 ALF (- Z)) P2) (* (+ J ALF) P3))
		    (1+ J)
		 )
	)
	(setq j (1+ j))
      )
      (setq pp (/ (- (* n p1) (* p2 (+ n alf))) z))
      (setq z1 z)
      (setq z (- z1 (/ p1 pp)))
      (setq its (1+ its))
    )
    (setq xi (cons (cons i z) xi))
    (setq yi (-	(/ (exp (- (Math::gammln (+ alf n)) (Math::gammln n)))
		   (* pp n p2)
		)
	     )
    )
    (setq wi (cons (cons i yi) wi))
    (setq i (1+ i))
  )
  (MATH::INT:Bind xi wi)
)

五、高斯-切比雪夫积分
此方法用于针对 sqrt(1-x^2)*f(x)型的积分有很高的效率
;;;=============================================================
;;; 高斯-切比雪夫积分                                           
;;; 此方法用于针对 sqrt(1-x^2)*f(x)型的积分有很高的效率。       
;;; 输入: foo 函数名,arg 除自变量x外的参数列表, n迭代次数。     
;;; 输出: 此类积分的数值.                                       
;;; 说明: 此积分法效率较高,n取值一般10次左右就达到有效浮点精度.
;;; http://mathworld.wolfram.com/Chebyshev-GaussQuadrature.html 
;;;=============================================================
(defun Math::INT:Gauss-Chebyshev (a b eps / FI FLAG I N S0 SX WI XI)
  (setq n 3)							;一般来说迭代6-5次左右就可以达到浮点计算精度
  (setq flag T)							;是否进行迭代
  (while (and (< n 1000) flag)					;迭代次数不超过100
    (setq wi (/ pi n))						;所有的权值均为pi/n
    (setq sx 0)
    (setq i 1)
    (repeat n
      (setq xi (cos (/ (* pi (- i 0.5)) n)))			;x 项
      (setq fi (MATH::INT:func xi))				;x 项的函数值
      (setq sx (+ sx fi))
      (setq i  (1+ i))
    )
    (setq sx (* wi sx))						;统一乘以权值
    (if (equal sx s0 eps)					;是否满足精度要求
      (setq flag nil)						;是则不再迭代
      (setq n (+ n n)						;否则迭代次数倍增
	    s0 sx						;存储积分值
      )
    )
  )
  sx
)


下面是一般函数的积分方法:
六 、龙贝格积分法:
;;;=============================================================
;;; 龙贝格积分法                                                
;;; Romberg Integration                                         
;;; 输入: 函数名--foo (用符号表示,一般形式是 (foo x a b c ...) 
;;;       参数表--args ,除去自变量x的其它的参数列表            
;;;       下区间--Ra                                            
;;;       上区间--Rb                                            
;;;       容许误差--eps                                         
;;; 输出: 所求函数在区间段的积分                                
;;;=============================================================
(defun MATH::INT:Romberg (a b eps / EP H I K M N P Q S X Y Y0)
  (setq h (- b a))
  (setq y nil)
  (setq i 0)
  (repeat 20
    (setq y (cons (cons i 0.0) y))
    (setq i (1+ i))
  )
  (setq y (reverse y))
  (setq y0 (* h (+ (MATH::INT:func a) (MATH::INT:func b)) 0.5))
  (setq y (cons (cons 0 y0) (cdr y)))
  (setq	m  1
	n  1
	ep (1+ eps)
  )
  (while (and (>= ep eps) (<= m 19))
    (setq p 0.0)
    (setq i 0)
    (repeat n
      (setq x (+ a (* (+ i 0.5) h)))
      (setq p (+ p (MATH::INT:func x)))
      (setq i (1+ i))
    )
    (setq p (/ (+ (cdar y) (* h p)) 2.0))
    (setq s 1.0)
    (setq k 1)
    (repeat m
      (setq s (+ s s s s))
      (setq q (/ (- (* s p) (cdr (assoc (1- k) y))) (1- s)))
      (setq y (subst (cons (1- k) p) (assoc (1- k) y) y))
      (setq p q)
      (setq k (1+ k))
    )
    (setq ep (abs (- q (cdr (assoc (1- m) y)))))
    (setq m (1+ m))
    (setq y (subst (cons (1- m) q) (assoc (1- m) y) y))
    (setq n (+ n n))
    (setq h (/ h 2.0))
  )
  q
)

七 辛普森积分法。
;;;=============================================================
;;; 辛普森积分法                                                
;;; Simpson Integration                                         
;;;=============================================================
(defun MATH::INT:Simpson (a b eps / EP H ITER K N P S1 S2 T1 T2 X)
  (setq n 1)
  (setq h (- b a))
  (setq t1 (* h (+ (MATH::INT:func a) (MATH::INT:func b)) 0.5))
  (setq s1 t1)
  (setq ep (1+ eps))
  (setq iter 0)
  (while (and (>= ep eps) (< iter 50))
    (setq p 0.0)
    (setq k 0)
    (repeat n
      (setq x (+ a (* (+ k 0.5) h)))
      (setq p (+ p (MATH::INT:func x)))
      (setq k (1+ k))
    )
    (setq t2 (/ (+ t1 (* h p)) 2.))
    (setq s2 (/ (- (* 4.0 t2) t1) 3.))
    (setq ep (abs (- s2 s1)))
    (setq t1 t2)
    (setq s1 s2)
    (setq n (+ n n))
    (setq h (/ h 2))
    (setq iter (1+ iter))
  )
  s2
)

八、变步长梯形求积分法 
;;;=============================================================
;;; 变步长梯形求积分法                                          
;;; Trapezoidal Integration 1                                   
;;;=============================================================
(defun MATH::INT:Trapezia (a b eps / H K N P S T1 T2 X iter)
  (setq n 1)
  (setq h (- b a))
  (setq t1 (* h (+ (MATH::INT:func a) (MATH::INT:func b)) 0.5))
  (setq p (1+ eps))
  (setq iter 0)
  (while (and (>= p eps) (< iter 100))
    (setq s 0)
    (setq k 0)
    (repeat n
      (setq x (+ a (* (+ k 0.5) h)))
      (setq s (+ s (MATH::INT:func x)))
      (setq k (1+ k))
    )
    (setq t2 (/ (+ t1 (* h s)) 2.))
    (setq p (abs (- t1 t2)))
    (setq t1 t2)
    (setq n (+ n n))
    (setq h (/ h 2))
    (setq iter (1+ iter))
  )
  t2
)


九、步长积分法。
;;;=============================================================
;;; 步长积分法                                                  
;;; Trapezoidal Integration 2                                   
;;;=============================================================
(defun MATH::INT:Trapzd (a b n / DEL IT SUM TNM X s)
  (if (= n 1)
    (setq s (* 0.5 (- b a) (+ (MATH::INT:func a) (MATH::INT:func b))))
    (progn
      (setq it 1)
      (repeat (- n 2)
	(setq it (lsh it 1))
      )
      (setq tnm it)
      (setq del (/ (- b a) tnm))
      (setq x (+ a (* 0.5 del)))
      (setq sum 0.0)
      (repeat it
	(setq sum (+ sum (MATH::INT:func x)))
	(setq x (+ x del))
      )
      (setq s (* 0.5 (+ s (/ (* (- b a) sum) tnm))))
    )
  )
)


十、自适应积分法
;;;=============================================================
;;; 自适应求积分法                                              
;;; Self-adapting Trapezia Integration                          
;;;=============================================================
(defun MATH::INT:Atrapezia (a b eps / F0 F1 H S T0 TT d)
  (setq d 1e-4)
  (setq h (- b a))
  (setq TT '(0. . 0.))
  (setq f0 (MATH::INT:func a))
  (setq f1 (MATH::INT:func b))
  (setq t0 (* h (+ f0 f1) 0.5))
  (car (MATH::INT:ppp a b h f0 f1 t0 eps d tt))
)
 
(defun MATH::INT:PPP (x0 x1 h f0 f1 t0 eps d tt / EPS1 F G P T1 T2 T3 X X2)
  (setq x (+ x0 (* h 0.5)))
  (setq f (MATH::INT:func x))
  (setq t1 (* h (+ f0 f) 0.25))
  (setq t2 (* h (+ f1 f) 0.25))
  (setq p (abs (- t0 t1 t2)))
  (if (or (< p eps) (< (* 0.5 h) d))
    (cons (+ (car tt) t1 t2) (cdr tt))
    (progn
      (setq g (* h 0.5))
      (setq eps1 (/ eps 1.4))
      (setq t3 (MATH::INT:ppp x0 x g f0 f t1 eps1 d tt))
      (setq t3 (MATH::INT:ppp x x1 g f f1 t2 eps1 d t3))
    )
  )
)

一些相关函数:
;;;=============================================================
;;; 功能: 第一类椭圆积分                                        
;;; 输入: Phi 和 k < 1                                          
;;; 输出: 所求椭圆积分值                                        
;;;=============================================================
(defun Math:Elliptic_Integral_1 (phi kCoff / )
  (defun MATH::INT:func (x / s)
    (setq s (* kCoff (sin x)))
    (/ 1.0 (sqrt (* (- 1 s) (1+ s))))
  )
  (MATH::INT:Romberg 0 phi 1e-15)
)
 
;;;=============================================================
;;; 功能: 第二类椭圆积分                                        
;;; 输入: Phi 和 k < 1                                          
;;; 输出: 所求椭圆积分值                                        
;;;=============================================================
(defun Math:Elliptic_Integral_2 (phi kCoff / )
  (defun MATH::INT:func (x /)
    (setq x (* kCoff (sin x)))
    (sqrt (* (- 1 x) (1+ x)))
  )
  (MATH::INT:Romberg 0 phi 1e-15)
)
 
;;;=============================================================
;;; 阶乘                                                        
;;;=============================================================
(defun Math::frac (n)
  (if (= n 0)
    1.
    (if	(= n 1)
      1.
      (* n (Math::frac (1- n)))
    )
  )
)
 
;;;=============================================================
;;; 伽玛函数                                                    
;;;=============================================================
(defun Math::gammln (xx / x y tmp ser j cof)
  (setq	cof '(76.18009172947146		 -86.50532032941677
	      24.01409824083091		 -1.231739572450155
	      0.1208650973866179e-2	 -0.5395239384953e-5
	     )
  )
  (setq x xx)
  (setq y x)
  (setq tmp (+ x 5.5))
  (setq tmp (- tmp (* (+ 0.5 x) (log tmp))))
  (setq ser 1.000000000190015)
  (setq j 0)
  (repeat 6
    (setq y (1+ y))
    (setq ser (+ ser (/ (nth j cof) y)))
    (setq j (1+ j))
  )
  (- (log (/ (* 2.5066282746310005 ser) x)) tmp)
)
 
;;;=============================================================
;;; 次函数用于比较排序                                          
;;;=============================================================
(defun MATH::INT:funcSort (e1 e2)
  (< (car e1) (car e2))
)
 
;;;=============================================================
;;; 组合两个系数                                                
;;;=============================================================
(defun MATH::INT:Bind (xi wi)
  (setq xi (vl-sort xi 'MATH::INT:funcSort))
  (setq wi (vl-sort wi 'MATH::INT:funcSort))
  (setq xi (mapcar 'cdr xi))
  (setq wi (mapcar 'cdr wi))
  (mapcar 'cons xi wi)
)



下面是用对话框创建的求积分法的例程:

;;;=============================================================
;;; 用各种方法求积分的程序                                      
;;;=============================================================
(defun C:Quadrature (/ ID OK DCL_FILE)
  (setq id (load_dialog (setq Dcl_File (MATH::INT:Write_Dcl))))	;从对话框中得到表达式
  (vl-file-delete Dcl_File)					;删除临时对话框文件
  (setq ok 2)
  (if (new_dialog "dcl_Integration" id)
    (progn
      (VL-CATCH-ALL-APPLY 'MATH::INT:GetSettings)		;读取默认数据
      (action_tile "help" "(MATH::INT:Help 1)")			;帮助
      (foreach k '(0 1 2 3 4 5 6 7 8)
	(setq k (strcat "K" (itoa k)))
        (action_tile k "(MATH::INT:OnBtn $key)")                ;按钮动作,对应相对的积分方法
      )		   
      (setq ok (start_dialog))
    )
  )
  (unload_dialog ID)
  (princ)
)
 
(defun C:JF (/)
  (VL-CATCH-ALL-APPLY 'C:Quadrature)
  (princ)
)
 
;;;=============================================================
;;; 从环境变量读取上次数据                                      
;;;=============================================================
(defun MATH::INT:GetSettings (/ data)
  (if (setq Data (getenv "Intergration"))
    (foreach k (read data)
      (set_tile (car k) (cdr k))
    )
  )
)
 
;;;=============================================================
;;; 检查对话框输入                                              
;;;=============================================================
(defun MATH::INT:CheckInput (symS symA symB symN / e f)
  (setq e (exp 1))
  (set symS (get_tile "F"))
  (set symA (MATH::INT:MyRead (get_tile "A")))
  (set symB (MATH::INT:MyRead (get_tile "B")))
  (set symN (MATH::INT:MyRead (get_tile "N")))
  (setq f (CAL:Expr2Func (eval s) 'MATH::INT:func '(x)))
  (apply 'and (mapcar 'eval '(F symS symA symB symN)))   
)
 
 
;;;=============================================================
;;; 按钮动作,对应相应的函数求积分                              
;;;=============================================================
(defun MATH::INT:OnBtn (key / DATA s N a b OldZIN m EPS RET tm0 map msg foo tmp idx)
  (setq m (VL-CATCH-ALL-APPLY 'MATH::INT:CheckInput '(s a b n)))
  (if (or (vl-catch-all-error-p m) (not m) (equal a b 1e-8))
    (if (vl-catch-all-error-p m) 
      (set_tile "info" (vl-catch-all-error-message m))
      (set_tile "info" "无效输入!")
    )
    (progn
      ;;如果精度过高,设置为15位的精度
      (if (> n 20)
	(setq n 15)
	(setq n (fix (abs n)))
      )
      ;;如果上区间小于下区间,则交换区间
      (if (< b a)
	(setq tmp a a b b tmp)
      )
      ;;记住对话框输入,用于下次
      (setq OldZIN (getvar "DIMZIN"))
      (setvar "DIMZIN" 8)
      (set_tile "N" (itoa n))
      (set_tile "A" (rtos a 2 20))
      (set_tile "B" (rtos b 2 20))
      (setq data (list (cons "F" s)
		       (cons "N" (itoa n))
		       (cons "A" (rtos a 2 20))
		       (cons "B" (rtos b 2 20))
		 )
      )
      (setvar "DIMZIN" OldZIN)
      (setenv "Intergration" (VL-PRIN1-TO-STRING data))
      ;;开始计算积分
      (setq eps (expt 0.1 n))
      (setq tm0 (getvar "TDUSRTIMER"))
      (setq map (MATH::INT:GetMethods))				;积分计算方法集
      (setq idx (atoi (substr key 2)))				
      (setq foo (nth idx map))					;获取计算积分的函数
      (setq ret (VL-CATCH-ALL-APPLY foo (list a b eps)))        ;获取积分值
      (if (vl-catch-all-error-p ret)
	(set_tile "info" (vl-catch-all-error-message ret))	;求解过程发生了错误
        (if (null ret)
	  (set_tile "info" "发生了错误,求值结果为空!")	
	  (progn
	    ;(MATH::INT:Bench 100 a b eps)
            (setq ret (rtos ret 2 20))
            (setq msg (get_attr key "label"))
            (setq msg (strcat msg "求的结果为:" ret))
            (set_tile (strcat "R" (itoa idx)) ret)             	;显示求解结果
            (set_tile "info" msg)				;显示求解结果
            (princ (strcat "\n" msg))				;打印求解结果
            (princ "\n费时:")
            (princ (* (- (getvar "TDUSRTIMER") tm0) 86400))
            (princ "秒.")
	  )
        )
      )
    )
  )
)
 
;;;=============================================================
;;; 各种积分测速                                                
;;;=============================================================
(defun MATH::INT:Bench (n a b eps)
  (UTI:BENCH
    n
    (list
      (list 'MATH::INT:Romberg a b eps)
      (list 'MATH::INT:Gauss-Legendre a b eps)
      (list 'MATH::INT:Simpson a b eps)
    )
  )
)
 
;;;=============================================================
;;; 积分方法集                                                  
;;;=============================================================
(defun MATH::INT:GetMethods ()
  '(MATH::INT:Romberg
    MATH::INT:Simpson
    MATH::INT:Atrapezia
    MATH::INT:Trapezia
    MATH::INT:Gauss-Legendre
    Math::INT:Gauss-Chebyshev
    Math::INT:Gauss-Laguerre
    Math::INT:Gauss-Hermite
    MATH::INT:Gauss-Jacobi
   )
)
 
;;;=============================================================
;;; 表达式求值,也可以用cal函数                                 
;;;=============================================================
(defun MATH::INT:MyRead (str / e)
  (setq e (exp 1))
  (CAL:Expr2Value str)
)
 
;;;=============================================================
;;; 帮助和说明: help and instruction                            
;;;=============================================================
(defun MATH::INT:Help (n)
  (if (= n 1)
    (if	(= "CHS" (getvar "Locale"))
      (alert
	"函数式只接受符号x为变量,不规范很可能出错!
	\n函数可以LISP内置的数学函数,也可以自定义函数!
	\n指数用^表示,+-*/表示加减乘除,乘号不能省略。
	\n程序能采用多种方法求积,一般来说龙贝格积分法最快。
	\n有什么问题email: highflybird@qq.com
	\n作者: highflybird 日期2019.07"
      )
      (alert
	"Standard expression only accepts \"x\" as a variale!
	\nThe fastest is Romberg Integration,the slowest is Trapezoidal rule(Be careful!).
	\nRecommendation:Don't set a high precision at first,promote it step by step.
	\nEspically for the Trapezoidal rule, It won't work well on some circumstances.
	\nIt's an Open Source Software. Thanks for your advice or bug reports.
	\nAuthor: highflybird  Email: highflybird@qq.com  Date:2019.07."
      )
    )
    (set_tile "info" "表达式非法或者空输入.")
  )
)

对话框的制作:
;;;=============================================================
;;; 输入对话框                                                  
;;;=============================================================
(defun UTI:Inputbox (/ str wcs ret)
  (setq	str "Function GetNumbers()
  	     GetNumbers=inputbox(\"请输入两个参数,中间用空格隔开:\",\"输入框\")
             End Function"
  )
  (if
    (or
      (setq wcs (vlax-create-object "Aec32BitAppServer.AecScriptControl.1"))
      (setq wcs (vlax-create-object "ScriptControl"))
    )
    (progn 
      (vlax-put-property wcs "language" "VBScript")
      (vlax-invoke wcs 'addcode str)
      (if (setq ret (vlax-invoke wcs 'run "GetNumbers"))
	(setq ret (strcat "(" ret ")")
	      ret (read ret)
	)
      )
      (vlax-release-object wcs)
      (if
	(and
	  (= 2 (length ret))
	  (or (= 'INT (type (car ret))) (= 'REAL (type (car ret))))
	  (or (= 'INT (type (cadr ret))) (= 'REAL (type (cadr ret))))
	)
	ret
      )
    )
  )
)
 
;;;=============================================================
;;; 写对话框到文件用于程序                                      
;;;=============================================================
(defun MATH::INT:Write_Dcl (/ Dcl_File file str)
  (setq Dcl_File (vl-filename-mktemp nil nil ".Dcl"))
  (setq file (open Dcl_File "w"))
  (princ
    "dcl_Integration : dialog {
	label = \"数值积分LISP版  v1.2\";
	: boxed_column {
          width = 60;
	  fixed_width = true;
	  : edit_box {
	    key=\"F\";
	    label= \"函数:\";
	  }
	  : row {
	    : edit_box {
	      key=\"A\";
	      label= \"下限:\";
	    }
	    : edit_box {
	      key=\"B\";
	      label= \"上限:\";
	    }      
	    : edit_box {
	      key=\"N\";
	      label = \"精确位数:\";
	      value = 8;
	      edit_width = 2;
	      fixed_width = true;
	    }
	  }
	  spacer_1;
	}
	: row {
	  : boxed_column {
	    label = \"计算方法:\";
	    : button {
	      key = \"K0\";
	      label = \"龙贝格积分法\";
	    }
	    : button {
	      key = \"K1\";
	      label = \"辛普森积分法\";
	    }
	    : button {
	      key = \"K2\";
	      label = \"自适应积分法\";
	    }
	    : button {
	      key = \"K3\";
	      label = \"变步长梯形法\";
	    }
	    : button {
	      key = \"K4\";
	      label = \"高斯-勒让德法\";
	    }
	    : button {
	      key = \"K5\";
	      label = \"高斯-切比雪夫法\";
	    }
	    : button {
	      key = \"K6\";
	      label = \"高斯-拉盖尔法\";
	    }
	    : button {
	      key = \"K7\";
	      label = \"高斯-埃尔米特法\";
	    }
	    : button {
	      key = \"K8\";
	      label = \"高斯-雅克比法\";
	    }
	    spacer;
	  }
	  : boxed_column {
	    width = 32;
	    fixed_width = true;
	    label = \"计算结果:\";
	    : text {
	      key = \"R0\";
	    }
	    : text {
	      key = \"R1\";
	    }
	    : text {
	      key = \"R2\";
	    }
	    : text {
	      key = \"R3\";
	    }
	    : text {
	      key = \"R4\";
	    }
	    : text {
	      key = \"R5\";
	    }
	    : text {
	      key = \"R6\";
	    }
	    : text {
	      key = \"R7\";
	    }
	    : text {
	      key = \"R8\";
	    }
	    spacer;
	  }
	}
	ok_cancel_help;
	//ok_cancel_help_errtile;
	: text {
	  key = \"info\";
	  label = \"Copyright \\u+00A9 2007-2019 Highflybird. All rights reserved.\";
	  width = 20;
	}
    }  "
    file
  )
  (close file)
  Dcl_File
)
 
;;;=============================================================
;;; 下面的就不用介绍了                                          
;;;=============================================================
(vl-load-com)
(if (= "CHS" (getvar "Locale"))
  (prompt "输入命令: JF")
  (prompt "Please enter: Quadrature")
)
(c:JF)
(princ)


一些应用,如下面网友的提问并解答:

http://bbs.mjtd.com/forum.php?mod=viewthread&tid=179809&extra=&highlight=%BB%FD%B7%D6&page=1

椭球体球缺的表面积计算
积分公式用LISP程序表达并计算出来
这是一个网上的
有结果和公式

用上面的介绍的一些函数可以得出其面积和体积:

;;;=============================================================
;;; 功能: 主程序,获取椭球冠的面积                              
;;; 输入: 无。                                                  
;;; 输出: 无。                                                  
;;;=============================================================
(defun c:test (/ a b c d h s0 s1)
  (initget 15)
  (setq a (getdist "\n输入椭球X半轴长: "))			
  (initget 15)
  (setq b (getdist "\n输入椭球y半轴长: "))
  (initget 15)
  (setq c (getdist "\n输入椭球z半轴长: "))
  (initget 15)
  (setq d (getdist "\n输入椭球冠的高度: "))			
  (if (<= d (+ c c))						;所求范围为小于椭球Z轴长
    (progn
      (setq h (/ (- c d) c))
      (setq s0 (ELL:GetSurfaceArea a b c 1))			;用龙贝格积分法计算其半个椭球面积
      (setq s2 (ELL:GetSurfaceArea a b c 1))			;用高斯-切比雪夫积分法计算其半个椭球面积
      (if (zerop h)						;为了防止被零除
	(setq s1 0 s2 0)
	(setq s1 (ELL:GetSurfaceArea a b c h) 			;用龙贝格积分法计算椭球带面积
	      s3 (ELL:GetSurfaceArea3 a b c h)			;用高斯-切比雪夫积分法计算椭球带面积
	)
      )
      (princ "\n用龙贝格法得到所求表面积是: ")
      (princ (rtos (- s0 s1) 2 20))				;两个面积相减,便是其椭球冠面积
 
      (princ "\n用高斯-切比雪夫积分法得到所求表面积是: ")
      (princ (rtos (- s2 s3) 2 20))				;两个面积相减,便是其椭球冠面积
 
      (princ "\n用龙贝格法所求的椭球的表面积是: ")
      (princ (rtos (+ s0 s0) 2 20))
 
      (princ "\n估算结果是: ")
      (princ (rtos (ELL:GetSurfaceArea2 a b c) 2 20))
 
      (princ "\n另外的结果是: ")
      (princ (rtos (ELL:GetSurfaceArea1 a b c) 2 20))
 
      (princ "\n高斯-切比雪夫积分的结果是: ")
      (princ (rtos (+ s2 s2) 2 20))
    )
  )
  (princ)
)
 
;;;=============================================================
;;; 功能: 获取椭球带的面积                                      
;;; 输入: 椭球的三个方向的轴径,和椭球带的高度                  
;;; 输出: 椭球带的面积                                          
;;;=============================================================
(defun ELL:GetSurfaceArea (a b c h / func)
  (defun Func (x a b c h / F G K S Y Z)
    (setq k (* a a))
    (setq s (sin x))
    (setq g (- 1 (* (- 1 (/ k b b)) s s)))
    (setq z (1- (/ k c c g)))
    (setq y (1+ (* z h h)))
    (if	(zerop z)
      (* (sqrt g) (1+ (sqrt y)))
      (if (< z 0)
	(setq z	(* h (sqrt (- z)))
	      f	(* (sqrt g) (+ (sqrt y) (/ (asin z) z)))
	)
	(setq z	(* h (sqrt z))
	      f	(* (sqrt g) (+ (sqrt y) (/ (asinh z) z)))
	)
      )
    )
  )
  (* 2 b c h (Math:Romberg 'func (list a b c h) 0 (* 0.5 pi) 1e-15))
)
 
;;;=============================================================
;;; 功能: 用椭圆函数计算椭球表面积                              
;;;=============================================================
(defun ELL:GetSurfaceArea1 (a b c / ll aa bb cc ac e1 e2 ph k f1 f2 s)
  (if (and (equal a b 1e-8) (equal b c 1e-8))
    (* 4 pi a)
    (setq ll (vl-sort (list a b c) '>)
	  a  (car ll)
	  b  (cadr ll)
	  c  (caddr ll)
	  aa (* a a)
	  bb (* b b)
	  cc (* c c)
	  ac (sqrt (- aa cc))
	  e1 (/ ac a)
	  e2 (/ (sqrt (- bb cc)) b)
	  ph (asin e1)
	  k  (/ e2 e1)
	  f1 (Math:Elliptic_Integral_1 ph k)
	  f2 (Math:Elliptic_Integral_2 ph k)
	  s  (* 2 pi (+ cc (/ (* b cc f1) ac) (* b ac f2)))
    )
  )
)
 
;;;=============================================================
;;; 功能: 估算椭球表面积                                        
;;;=============================================================
(defun ELL:GetSurfaceArea2 (a b c)
  (* 4
     pi
     (expt
       (/ (+ (expt (* a b) 1.6075)
	     (expt (* b c) 1.6075)
	     (expt (* c a) 1.6075)
	  )
	  3
       )
       (/ 1 1.6075)
     )
  )
)
 
 
 
;;;=============================================================
;;; 几个相关数学函数                                            
;;;=============================================================
(defun asin (x) (atan x (sqrt (- 1 (* x x))))) 			;反正弦函数
(defun asinh (x) (log (+ x (sqrt (1+ (* x x))))))		;反双曲正弦函数=log(x+sqrt(x*x+1))
 
(vl-load-com)
(princ "\n运行命令是:Test")
(princ)


至于一些积分函数有什么区别,我就不在这里一一介绍了。

希望这些代码对你有帮助。引用请注明来源。

LISP的随机函数的实现

对Autolisp来说,没有生成随机数的函数,一般来说,有如下方法可以产生它:

方法一:用同余函数实现,这是纯LISP的方法。

;;;rand function-some code from Lee Mac,--thanks.
(defun LM:random (Var / x)
  (/ (setq x	4294967296.0
	   seed	(rem (1+ (* 1664525.0
			    (cond (seed)
				  ((getvar var))
			    )
			 )
		     )
		     x
		)
     )
     x
  )
)
(defun Rand (a b var)
  (fix (+ a (* (LM:random var) (- b a -1))))
)

方法二:用系统的相关的变量如cputicks,TDUSRTIMER等实现:

此方法可能会产生一些不是很随机的数。

(defun c:random	()
  (setq seed (getvar "TDUSRTIMER"))
  (setq seed (- seed (fix seed)))
  (rem (* seed 86400) 1)
)

方法三:利用windows系统com中产生随机数。

当前实现Random类基于 Donald E.Knuth 删减随机数生成器算法的修改版本。 有关详细信息,请参阅 D.e。 Knuth。 计算机编程,卷 2 的技术:Seminumerical 算法。 Addison-wesley,Reading,MA,第三个版本,1997年。

(setq rndObj (vlax-create-object "System.Random"))
  (setq x (vlax-invoke rndobj 'nextdouble))
  (setq y (vlax-invoke rndobj 'next))
 
  (defun sysRnd (x y)
    (+ x (* (vlax-invoke rndobj 'nextdouble) (- y x)))
  )

方法四:用vbscript 或者Jscirpt等内部的随机函数产生,然后在lisp用调用脚本语言。

(if
    (or
      (setq wcs (vlax-create-object   "Aec32BitAppServer.AecScriptControl.1"))
      (setq wcs (vlax-create-object "ScriptControl"))
    )
    (vlax-put-property wcs "language" "VBScript")
  )
  (setq str "Randomize
             Function Rand(x,y)\n
                 Rand=x+Rnd*(y-x)\n
             End Function"
  )        
  (vlax-invoke wcs 'addcode str)
 
  (defun wcsRnd(a b)
    (vlax-invoke wcs 'run "rand" a b)
  )

上面的代码经过测试验证,用同余的速度最快,vbs最慢,system.random中间,但同余法的效率并未超过system.random很多,仅仅是几倍的差距。而system.random虽也是伪随机的,但是正确性比同余法高。

 

 



一元二次、三次、四次方程求解和复数的运算

对一元三次或者四次方程,是有数学公式求精确解的,可以不用迭代法。参考了维基的上的方法,现在我贴出一元二次、三次或者四次方程的LISP求解方法。使得在求解效率可以得到极大提高。
注明: 因为这几个方程的解有可能是复数,所以我对每个解都用表的形式来列出。
如果这个表的第二项为0,那么这个解是实数,否则是复数。
譬如 :{x^4+3x^3+7x^2+2x-5 = 0}

(Math:Quartic_Equation 1 3 7 2 -5)

==》((-1.19281 -2.21406) (-1.19281 2.21406) (-1.24789 0) (0.633498 0))

意味着这个方程有两个实数解:-1.24789 , 0.633498
两个虚数解:-1.19281-2.21406 i ,-1.19281+2.21406i
另外在末尾附上验算测试函数。

解一元二次方程的LISP源码:

  1. ;;;=============================================================
  2. ;;;一元二次方程的解                                     	
  3. ;;;f(x) = a*x^2+b*x+c = 0                               	
  4. ;;;Input: the coefficients a, b, c are real numbers     	
  5. ;;;Output: when a /= 0,one or Two solutions,all of them like 	
  6. ;;;        this: (x y),means: x + y * i,if it's a real number,  
  7. ;;;        then y = 0.Otherwise ,  return a real number or nil  
  8. ;;;Ref: http://en.wikipedia.org/wiki/Quadratic_equation 	
  9. ;;;=============================================================
  10. (defun Math:Quadratic_Equation (a b c / d e g)
  11.   (if (zerop a)
  12.     (if (not (zerop b))
  13.       (list (list (/ (- c) (float b)) 0))
  14.     )
  15.     (progn
  16.       (setq a (float a))
  17.       (if (not (equal a 1 1e-14))
  18.         (setq b (/ b a)
  19.               c (/ c a)
  20.         )
  21.       )
  22.       (setq d (- (* b b) (* 4 c)))
  23.       (setq e (* b -0.5))
  24.       (cond
  25.         ( (equal d 0 1e-14)
  26.           (list (list e 0) (list e 0))
  27.         )
  28.         ( (> d 0)
  29.           (setq g (* (sqrt d) -0.5))
  30.           (list (list (- e g) 0) (list (+ e g) 0))
  31.         )
  32.         ( (< d 0)
  33.           (setq g (* (sqrt (- d)) -0.5))
  34.           (list (list e (- g)) (list e g))
  35.         )
  36.       )
  37.     )
  38.   )
  39. )

解一元三次方程的LISP源码:

  1. ;;;=============================================================
  2. ;;;一元三次方程的解                                     	
  3. ;;;f(x) = a*x^3+b*x^2+c*x+d = 0                         	
  4. ;;;Input: the coefficients a, b, c, d are real numbers  	
  5. ;;;Output: when a /= 0,Three solutions,all of them like this: 	
  6. ;;;        (x y),means: x+y*i,if it's a real number,then y = 0; 
  7. ;;;        otherwise goto "Math:Quadratic_Equation"     	
  8. ;;;Ref: http://en.wikipedia.org/wiki/Cubic_function     	
  9. ;;;=============================================================
  10. (defun Math:Cubic_Equation (a b c d / e f g h u w x y)
  11.   (cond
  12.     ( (zerop a)
  13.       (Math:Quadratic_Equation b c d)
  14.     )
  15.     ( (zerop d)
  16.       (cons '(0 0) (Math:Quadratic_Equation a b c))
  17.     )
  18.     (t
  19.       (setq b (/ b a 3.0))
  20.       (setq c (/ c a 6.0))
  21.       (setq d (/ d a 2.0))
  22.  
  23.       (setq e (- (* b (- (+ c c c) (* b b))) d))        	;Alpha
  24.       (setq f (- (* b b) c c))                          	;Beta
  25.       (setq g (- (* e e) (* f f f)))                    	;Delta,The nature of the roots
  26.       (cond
  27.         ( (equal g 0 1e-14)
  28.           (setq u (MATH:CubeRoot e))
  29.           (list (list (- (+ u u) b) 0.0)
  30.                 (list (- (+ b u)) 0.0)
  31.                 (list (- (+ b u)) 0.0)
  32.           )
  33.         )
  34.         ( (> g 0)
  35.           (setq h (sqrt g))
  36.           (setq u (MATH:CubeRoot (+ e h)))
  37.           (setq w (MATH:CubeRoot (- e h)))
  38.           (setq x (- (+ b (* (+ u w) 0.5))))
  39.           (setq y (* (sqrt 3) (- u w) 0.5))
  40.           (list (list (+ u w (- b)) 0)
  41.                 (list x y)
  42.                 (list x (- y))
  43.           )
  44.         )
  45.         ( (< g 0)
  46.           (setq x (/ e f (sqrt f)))
  47.           (setq y (sqrt (abs (- 1 (* x x)))))
  48.  
  49.           (setq u (/ (atan y x) 3))
  50.           (setq w (/ (+ pi pi) 3))
  51.           (setq h (* 2 (sqrt f)))
  52.           (list (list (- (* (cos u) h) b) 0)
  53.                 (list (- (* (cos (+ u w)) h) b) 0)
  54.                 (list (- (* (cos (- u w)) h) b) 0)
  55.           )
  56.         )
  57.       )
  58.     )
  59.   )
  60. )

解一元四次方程的解:

  1. ;;;=============================================================
  2. ;;;一元四次方程的解                                     	
  3. ;;;f(x) = a*x^4+b*x^3+c*x^2+d*x+e= 0                    	
  4. ;;;Input: the coefficients a,b,c,d,e are real numbers.  	
  5. ;;;Output: when a /= 0,Three solutions,all of them like this: 	
  6. ;;;        (x y),means: x+y*i,if it's a real number,then y = 0, 
  7. ;;;        otherwise goto "Math:Quadratic_Equation".            
  8. ;;;Ref: http://en.wikipedia.org/wiki/Quartic_function   	
  9. ;;;=============================================================
  10. (defun Math:Quartic_Equation (a b c d e / B2 B3 B4 EPS F G H P Q R V W X Y Z S)
  11.   (setq eps 1e-8)
  12.   (cond
  13.     ( (equal a 0 eps)
  14.       (Math:Cubic_Equation b c d e)
  15.     )
  16.     ( (equal e 0 eps)
  17.       (cons '(0 0) (Math:Cubic_Equation a b c d))
  18.     )
  19.     ( (and (equal b 0 eps) (equal d 0 eps))
  20.       (foreach x (Math:Quadratic_Equation a c e)
  21.         (foreach y (Math:CSqrt x)
  22.           (setq s (cons y s))
  23.         )
  24.       )
  25.     )  
  26.     ( (setq a (float a))
  27.       (setq b (/ b a))
  28.       (setq c (/ c a))
  29.       (setq d (/ d a))
  30.       (setq e (/ e a))
  31.       (setq b2 (* b b))
  32.       (setq b3 (* b2 b))
  33.       (setq b4 (* b3 b))
  34.  
  35.       (setq f (+ (* b2 -0.375) c))                      	;alpha
  36.       (setq g (+ (* b3 0.125) (* b c -0.5) d))          	;beta
  37.       (setq h (+ (* -0.01171875 b4) (* 0.0625 c b2) (* -0.25 b d) e))
  38.  
  39.       (if (equal g 0 eps)
  40.         (progn 
  41.           (setq x (* b -0.25))
  42.           (mapcar
  43.             (function (lambda (e) (list (+ (car e) x) (cadr e))))
  44.             (Math:Quartic_Equation 1 0 f 0 h)           	;Recursion 
  45.           )
  46.         )
  47.         (progn 
  48.           (setq p (- (* f f 2) h))
  49.           (setq q (- (* f f f 0.5) (* f h 0.5) (* g g 0.125)))
  50.           (setq r (Math:Cubic_Equation 1 (* 2.5 f) p q))
  51.           (while (not (equal (cadar r) 0 eps))
  52.             (setq r (cdr r))
  53.           )
  54.           (setq r (caar r))
  55.           (setq w (sqrt (abs (+ f r r))))
  56.           (foreach i (list + -)
  57.             (foreach j (list + -)
  58.               (setq x (i (* b -0.25) (* w 0.5)))
  59.               (setq y (- (+ f f f r r (i (/ g 0.5 w)))))
  60.               (if (< y 0)
  61.                 (setq z (list x (j (* (sqrt (- y)) 0.5))))
  62.                 (setq z (list (+ x (j (* (sqrt y) 0.5))) 0)) 
  63.               )
  64.               (setq S (cons z S))
  65.             )
  66.           )
  67.         )
  68.       )
  69.     )
  70.   )
  71. )

复数相关代码:

  1. ;;;=============================================================
  2. ;;;立方根(为解三次和四次方程编写,因expt函数第一个参数不能为负) 
  3. ;;;输入:x 实数                                                 
  4. ;;;输出:x 的立方根                                             
  5. ;;;=============================================================
  6. (defun MATH:CubeRoot (x)
  7.   (if (< x 0)
  8.     (- (expt (- x) 0.33333333333333333333333))
  9.     (expt x 0.33333333333333333333333)
  10.   )
  11. )
  12.  
  13. ;;;=============================================================
  14. ;;;复数加复数                                                   
  15. ;;;输入:Z1--复数,Z2--复数                                     
  16. ;;;输出:复数跟实数相加的结果                                   
  17. ;;;=============================================================
  18. (defun Math:C+C (z1 z2)
  19.   (mapcar '+ z1 z2)
  20. )
  21.  
  22. ;;;=============================================================
  23. ;;;复数减复数                                                   
  24. ;;;输入:Z1--复数,Z2--复数                                     
  25. ;;;输出:复数跟实数相减的结果                                   
  26. ;;;=============================================================
  27. (defun math:C-C (z1 z2)
  28.   (mapcar '- Z1 Z2)
  29. )
  30.  
  31. ;;;=============================================================
  32. ;;;复数乘复数                                                   
  33. ;;;输入:Z1--复数,Z2--复数                             	
  34. ;;;输出:复数跟实数相乘的结果                           	
  35. ;;;=============================================================
  36. (defun Math:C*C (Z1 Z2)
  37.   (list
  38.     (- (* (car Z1) (car z2)) (* (cadr z1) (cadr z2)))
  39.     (+ (* (car Z1) (cadr Z2)) (* (cadr z1) (car Z2)))
  40.   )
  41. )
  42.  
  43. ;;;=============================================================
  44. ;;;复数除以复数                                         	
  45. ;;;输入:Z1--复数,Z2--复数                             	
  46. ;;;输出:复数跟实数相除的结果                           	
  47. ;;;=============================================================
  48. (defun Math:C/C (Z1 Z2 / a b c d N)
  49.   (setq a (car  z1))
  50.   (setq b (cadr z1))
  51.   (setq c (car  z2))
  52.   (setq d (cadr z2))
  53.   (setq N (float (+ (* c c) (* d d))))
  54.   (if (not (zerop N))
  55.     (list
  56.       (/ (+ (* a c) (* b d)) N)
  57.       (/ (- (* b c) (* a d)) N)
  58.     )
  59.   )
  60. )
  61.  
  62. ;;;=============================================================
  63. ;;;复数加实数                                           	
  64. ;;;输入:Z --复数,R --实数                             	
  65. ;;;输出:复数跟实数相加的结果                           	
  66. ;;;=============================================================
  67. (defun Math:C+R (Z R)
  68.   (list (+ (car z) R) (cadr z))
  69. )
  70.  
  71. ;;;=============================================================
  72. ;;;复数减实数                                           	
  73. ;;;输入:Z --复数,R --实数                             	
  74. ;;;输出:复数跟实数相减的结果                           	
  75. ;;;=============================================================
  76. (defun Math:C-R (Z R)
  77.   (list (- (car z) R) (cadr z))
  78. )
  79.  
  80. ;;;=============================================================
  81. ;;;复数乘实数                                           	
  82. ;;;输入:Z --复数,R --实数                             	
  83. ;;;输出:复数跟实数相乘的结果                           	
  84. ;;;=============================================================
  85. (defun Math:C*R (Z R)
  86.   (list (* (car z) R) (* (cadr z) R))
  87. )
  88.  
  89. ;;;=============================================================
  90. ;;;复数除以实数                                         	
  91. ;;;输入:Z --复数,R --实数                             	
  92. ;;;输出:复数跟实数相除的结果                           	
  93. ;;;=============================================================
  94. (defun Math:C/R (Z R)
  95.   (if (not (zerop R))
  96.     (list (/ (car z) R) (/ (cadr z) R))
  97.   )
  98. )
  99.  
  100. ;;;=============================================================
  101. ;;;复数的模                                             	
  102. ;;;输入:Z --复数                                       	
  103. ;;;输出:复数的模,即复数代表的矢量长度                 	
  104. ;;;=============================================================
  105. (defun MATH:CNormal (Z)
  106.   (distance '(0 0) Z)
  107. )
  108.  
  109. ;;;=============================================================
  110. ;;;复数的平方根                                         	
  111. ;;;输入:Z --复数                                       	
  112. ;;;输出:复数的平方根,以复数形式表达,有两个解         	
  113. ;;;=============================================================
  114. (defun Math:CSqrt (Z / x y a r)
  115.   (setq x (car z))
  116.   (setq y (cadr z))
  117.   (if (equal y 0 1e-14)
  118.     (if (> x 0)
  119.       (list (list (setq x (sqrt x)) 0) (list (- x) 0))
  120.       (list (list 0 (setq y (sqrt (- x)))) (list 0 (- y)))
  121.     )
  122.     (progn
  123.       (setq a (* (atan y x) 0.5))
  124.       (setq r (sqrt (distance '(0 0) Z)))
  125.       (setq x (* r (cos a)))
  126.       (setq y (* r (sin a)))
  127.       (list (list x y) (list (- x) (- y)))
  128.     )
  129.   )
  130. )
  131.  
  132. ;;;=============================================================
  133. ;;;复数的方根                                           	
  134. ;;;输入:Z --复数, n --正整数                           	
  135. ;;;输出:复数的n次方根,以复数形式表达,有n个解。       	
  136. ;;;=============================================================
  137. (defun MATH:CRoot (Z n / r a b i s 2Pi)
  138.   (if (and (= (type n) 'INT) (> n 0))
  139.     (progn
  140.       (setq r (expt (distance '(0 0) z) (/ 1. n)))
  141.       (setq a (atan (cadr z) (car z)))
  142.       (setq i 0)
  143.       (setq 2Pi (+ pi pi))
  144.       (repeat n
  145.         (setq b (/ (+ a (* i 2Pi)) n))
  146.         (setq s (cons (list (* r (cos b)) (* r (sin b))) s))
  147.         (setq i (1+ i))
  148.       )
  149.       (reverse s)
  150.     )
  151.   )
  152. )
  153.  
  154. ;;;=============================================================
  155. ;;;复数的实幂指数                                       	
  156. ;;;输入:Z --复数, x -- 实数                            	
  157. ;;;输出:复数Z 的x次幂,以复数形式表达。                	
  158. ;;;=============================================================
  159. (defun MATH:CRPower (Z x / a r)
  160.   (setq a (atan (cadr Z) (car Z)))
  161.   (setq r (expt (distance '(0 0) Z) x))
  162.   (list (* r (cos (* a x))) (* r (sin (* a x))))
  163. )
  164.  
  165. ;;;=============================================================
  166. ;;;复数的复幂指数                                       	
  167. ;;;输入:Z1 --复数, Z2 -- 复数                          	
  168. ;;;输出:复数Z1的Z2次幂,以复数形式表达。               	
  169. ;;;=============================================================
  170. (defun MATH:CZPower (Z1 Z2 / a r x y k)
  171.   (if (> (setq r (distance '(0 0) Z1)) 1e-20)
  172.     (progn 
  173.       (setq a (atan (cadr Z1) (car Z1)))
  174.       (setq x (car z2))
  175.       (setq y (cadr z2))
  176.       (setq k (log r))
  177.       (setq r (exp (- (* k x) (* y a))))
  178.       (setq a (+ (* k y) (* x a)))
  179.       (list (* r (cos a)) (* r (sin a)))
  180.     )
  181.   )
  182. )
  183.  
  184. ;;;=============================================================
  185. ;;;复数的自然对数                                        	
  186. ;;;输入:Z --复数                                       	
  187. ;;;输出:复数Z 的自然对数,以复数形式表达。             	
  188. ;;;=============================================================
  189. (defun MATH:CExp (Z / r)
  190.   (if (> (setq r (distance '(0 0) Z)) 1e-20)
  191.     (list (log r) (atan (cadr Z) (car Z)))
  192.   )
  193. )
  194.  
  195. ;;;=============================================================
  196. ;;;复数的余弦                                           	
  197. ;;;输入:Z --复数                                          	
  198. ;;;输出:复数Z 的余弦,以复数形式表达。                 	
  199. ;;;=============================================================
  200. (defun MATH:CCOS (Z / x y c s u v)
  201.   (setq x (car z))
  202.   (setq y (cadr z))
  203.   (if (equal y 0 1e-20)
  204.     (list (cos x) 0)
  205.     (progn
  206.       (setq c (* 0.5 (cos x)))
  207.       (setq s (* 0.5 (sin x)))
  208.       (setq u (exp y))
  209.       (setq v (exp (- y)))
  210.       (list (* c (+ v u)) (* s (- v u)))
  211.     )
  212.   )
  213. )
  214.  
  215. ;;;=============================================================
  216. ;;;复数的正弦                                          		
  217. ;;;输入:Z --复数                                       	
  218. ;;;输出:复数Z 的正弦,以复数形式表达。                 	
  219. ;;;=============================================================
  220. (defun MATH:CSIN (Z / x y c s u v)
  221.   (setq x (car z))
  222.   (setq y (cadr z))
  223.   (if (equal y 0 1e-20)
  224.     (list (sin x) 0)
  225.     (progn
  226.       (setq s (* 0.5 (sin x)))
  227.       (setq c (* 0.5 (cos x)))
  228.       (setq u (exp (- y)))
  229.       (setq v (exp y))
  230.       (list (* s (+ v u)) (* c (- v u)))
  231.     )
  232.   )
  233. )
  234.  
  235. ;;;=============================================================
  236. ;;;复数的正切                                           	
  237. ;;;输入:Z --复数                                       	
  238. ;;;输出:复数Z 的正切,以复数形式表达。                 	
  239. ;;;=============================================================
  240. (defun MATH:CTAN (Z / x y c s u v)
  241.   (MATH:C/C (MATH:CSIN Z) (MATH:CCOS Z))
  242. )
  243.  
  244. ;;;=============================================================
  245. ;;;实系数的复数多项式运算                               	
  246. ;;;输入:Z --复数, RealNumbers --实系数列表             	
  247. ;;;输出:根据实系数多项式复数方程值,用复数表示。       	
  248. ;;;=============================================================
  249. (defun MATH:CReal_Polynomial (Z RealNumbers / f)
  250.   (cond
  251.     ( (cdr RealNumbers) 
  252.       (setq f (MATH:C*R Z (car RealNumbers)))
  253.       (repeat (- (length RealNumbers) 2)
  254.         (setq RealNumbers (cdr RealNumbers))
  255.         (setq f (MATH:C+R f (car RealNumbers)))
  256.         (setq f (MATH:C*C f Z))
  257.       )
  258.       (setq f (MATH:C+R f (cadr RealNumbers)))
  259.     )
  260.     ( (car RealNumbers) (list (car RealNumbers) 0))
  261.   )
  262. )
  263.  
  264. ;;;=============================================================
  265. ;;;虚系数的复数多项式运算                               	
  266. ;;;输入:Z --复数, ImaginaryNumbers--复系数列表         	
  267. ;;;输出:根据复系数多项式复数方程值,用复数表示。       	
  268. ;;;=============================================================
  269. (defun MATH:CImaginary_Polynomial (Z ImaginaryNumbers / f)
  270.   (if (setq f (car ImaginaryNumbers))
  271.     (progn
  272.       (repeat (1- (length ImaginaryNumbers))
  273.         (setq ImaginaryNumbers (cdr ImaginaryNumbers))
  274.         (setq f (MATH:C*C f Z)) 
  275.         (setq f (MATH:C+C f (car ImaginaryNumbers)))
  276.       )
  277.       f
  278.     )
  279.   )
  280. )

测试代码:

  1. ;;;=============================================================
  2. ;;;以下程序为测试用:                                    	
  3. ;;;Test for "Math:Quartic_Equation"                     	
  4. ;;;If the difference between the Coefficients is big,it will get
  5. ;;;some float point error.                   			
  6. ;;;=============================================================
  7. (defun C:t4 (/ a b c d e S z z1 z2 z3 z4)
  8.   (initget 1)
  9.   (setq a (getreal "\nCoefficient a:"))
  10.   (initget 1)
  11.   (setq b (getreal "\nCoefficient b:"))
  12.   (initget 1)
  13.   (setq c (getreal "\nCoefficient c:"))
  14.   (initget 1)
  15.   (setq d (getreal "\nCoefficient d:"))
  16.   (initget 1)
  17.   (setq e (getreal "\nCoefficient e:"))
  18.   ;(MISC:test 1000 '((Math:Quartic_Equation a b c d e)))
  19.   (setq S (Math:Quartic_Equation a b c d e))            	;get the solutions
  20.   (foreach z S
  21.     (princ "\n")
  22.     (princ (mapcar '(lambda (x) (rtos x 2 20)) z))      	;print them.
  23.   )
  24.   (foreach z S                                          	;Check every solution
  25.     (setq f (MATH:CReal_Polynomial z (list a b c d e)))
  26.     (if (not (equal f '(0 0) 1e-6))                     	;if f(z) /= 0.0,maybe it's caused by floating point operation;
  27.       (alert
  28.         (strcat
  29.           "Maybe it's a Wrong solution: f("
  30.           (vl-princ-to-string z)
  31.           ") = "
  32.           (VL-PRINC-TO-STRING f)
  33.         )
  34.       )
  35.     )
  36.   )
  37.   (princ)
  38. )
  39.  
  40. (defun C:t3 (/ a b c d e S z z1 z2 z3 z4)
  41.   (initget 1)
  42.   (setq a (getreal "\nCoefficient a:"))
  43.   (initget 1)
  44.   (setq b (getreal "\nCoefficient b:"))
  45.   (initget 1)
  46.   (setq c (getreal "\nCoefficient c:"))
  47.   (initget 1)
  48.   (setq d (getreal "\nCoefficient d:"))
  49.   (UTI:BENCH
  50.     20000 
  51.     (list
  52.       (list 'Math:Cubic_Equation a b c d)
  53.       (list 'Math:Cubic_Equation1 a b c d)
  54.     )
  55.   )
  56.   (foreach n '(Math:Cubic_Equation Math:Cubic_Equation1)
  57.     (setq S (apply n (list a b c d)))            		;get the solutions
  58.     (foreach z S
  59.       (princ "\n")
  60.       (princ (mapcar '(lambda (x) (rtos x 2 20)) z))      	;print them.
  61.     )
  62.     (foreach z S                                          	;Check every solution
  63.       (setq f (MATH:CReal_Polynomial z (list a b c d)))
  64.       (if (not (equal f '(0 0) 1e-6))                     	;if f(z) /= 0.0,maybe it's caused by floating point operation;
  65.         (alert
  66.           (strcat
  67. 	    "Use "
  68. 	    (VL-PRINC-TO-STRING n)
  69.             ",Maybe it's a Wrong solution: f("
  70.             (vl-princ-to-string z)
  71.             ") = "
  72.             (VL-PRINC-TO-STRING f)
  73.           )
  74.         )
  75.       )
  76.     )
  77.   )
  78.   (princ)
  79. )
  80. (princ "\nThe command is : test.")
  81. (vl-load-com)
  82. (princ)

曲线的转弯半径和曲率

在下面的这个帖子中讨论了椭圆的曲率和转弯半径
http://bbs.mjtd.com/thread-62980-1-1.html
现在我把这个主题深化一下,讨论一下曲线的两个函数:
vlax-curve-getSecondDeriv
vlax-curve-getFirstDeriv
这两个函数是什么意思呢?
我们考察AutoCAD里面的曲线类,主要是圆,椭圆,弧和样条曲线,多段线由这几种组合而成。
椭圆和样条曲线实际上都是由参数形成,因此,对于这类曲线,它们每点的坐标可以由参数方程表达:
譬如椭圆 x=a*cos(t); y=b*sin(t);
样条曲线也有方程,假设样条曲线的参数方程为: X= f(t);
Y=g(t);
因此可以对参数方程求导,得到每一点的切线矢量,曲线上每一点对应于一个参数 t0 ,
这个切线矢量的 的X值就是 f(t)在t0处的一阶导数,Y值就是g(t)在t0处的一阶导数,
即( f'(t0), g'(t0),0)
因而我们就理解了vlax-curve-getFirstDeriv 函数返回值的意义,
对于vlax-curve-getSecondDeriv的意义类似,只不过这次换成了二阶导数。
即( f”(t0), g”(t0),0)

那么如何求样条曲线或者椭圆的每一点的曲率及其半径呢?
根据参数方程的求曲率公式:

 

下面是求值代码:

;;;=============================================================
;;;功能: 获取椭圆上一点处的曲率和离心半径(适用于单次利用此法)
;;;参数: 椭圆实体和椭圆上的一点
;;;返回: 此处离心圆圆心、离心半径及其曲率(离心率)
;;;说明: 如果要在CAD中几何作图,可以参考此贴:
;;;      http://bbs.mjtd.com/thread-62980-1-1.html
;;;=============================================================
(defun ELL:GetCurvature (en pt / obj a b px x y v1 v2 rad cen)
  (setq obj (vlax-ename->vla-object en))
  (setq a   (vla-get-MajorRadius obj))                          ;椭圆的半长轴
  (setq b   (vla-get-MinorRadius obj))				;椭圆的半短轴
  (setq pt  (vlax-curve-getclosestpointto en pt))		;保证此点在椭圆上
  (setq px  (vlax-curve-getParamAtPoint en pt))			;此点的椭圆参数
  (setq v1  (vlax-curve-getFirstDeriv en px))			;此点的一阶矢量
  (setq v2  (list (- (cadr v1)) (car v1) (caddr v1)))		;此点的切线矢量
  (setq x   (* a (sin px)))
  (setq y   (* b (cos px)))
  (setq rad (/ (expt (+ (* x x) (* y y)) 1.5) (* a b)))		;得到转弯半径
  (setq cen (polar pt (angle '(0 0 0) v2) rad))			;圆心
  (list cen rad)						;圆心及半径
)
;;;=============================================================
;;; 利用参数方程的求椭圆的曲率及其离心半径(适用于多次利用此法)
;;; 功能: 获取椭圆上一点处的转弯半径和离心圆圆心
;;; 参数: 椭圆实体和曲线上的一点
;;; 返回: 此处离心圆圆心、离心半径
;;;=============================================================
(defun ELL:GetCurvature1 (en a b pt / px v1 v2 x y rad cen)
  (setq pt  (vlax-curve-getclosestpointto en pt))		;保证此点在椭圆上
  (setq px  (vlax-curve-getParamAtPoint en pt))			;此点的椭圆参数
  (setq v1  (vlax-curve-getFirstDeriv en px))			;此点的一阶矢量
  (setq v2  (list (- (cadr v1)) (car v1) (caddr v1)))		;此点的切线矢量
  (setq x   (* a (sin px)))
  (setq y   (* b (cos px)))
  (setq rad (/ (expt (+ (* x x) (* y y)) 1.5) (* a b)))		;得到转弯半径
  (setq cen (polar pt (angle '(0 0 0) v2) rad))			;圆心
  (list cen rad)						;圆心及半径
)
;;;=============================================================
;;; 一般平面曲线参数方程的曲率离心公式
;;; 功能: 获取曲线上一点处的离心半径和离心圆圆心
;;; 参数: 曲线实体和曲线上的一点
;;; 返回: 此处离心圆圆心、离心半径
;;;=============================================================
(defun CUR:GetCurvature (en pt / ob px v1 v2 v3 x1 y1 x2 y2 cen rad d1 d2)
  (setq ob (vlax-ename->vla-object en))
  (setq pt (vlax-curve-getclosestpointto en pt))                ;保证此点在曲线上
  (setq px (vlax-curve-getParamAtPoint en pt))                  ;此点的曲线参数
  (setq v1 (vlax-curve-getFirstDeriv en px))                    ;此点的一阶矢量
  (setq v2 (vlax-curve-getSecondDeriv en px))			;此点的二阶矢量
  (setq v3 (list (- (cadr v1)) (car v1) (caddr v1)))            ;此点的切线矢量
  (setq x1 (car  v1))                                           ;一阶导数的 X值
  (setq y1 (cadr v1))						;一阶导数的 Y值
  (setq x2 (car  v2))						;二阶导数的 X值
  (setq y2 (cadr v2))						;二阶导数的 Y值
  (setq d1 (expt (+ (* y1 y1) (* x1 x1)) 1.5))
  (setq d2 (- (* x1 y2) (* x2 y1)))				;转弯内外的判定
  (if (/= d2 0)                                                 ;如果不为直线段
    (progn
      (setq rad (/ d1 d2))
      (if (vlax-method-applicable-p ob 'GetBulge)               ;如果为多段线(含圆弧)
	(if (< (vla-GetBulge ob (fix px)) 0)			;如果此段凸度小于0
	  (setq rad (- rad))
        )
      )
      (list (polar pt (angle '(0 0 0) v3) rad) (abs rad))	;圆心及半径
    )
  )
)

一些测试代码:

;;;=============================================================
;;;测试程序1: 获取曲线一点处的曲率和离心半径
;;;=============================================================
(defun c:tt1 (/ ent pnt ret)
  (setq ent (car (entsel "\n请选取曲线:")))
  (while (setq pnt (getpoint "\n点取一点"))
    (setq pnt (trans pnt 1 0))
    (setq pnt (vlax-curve-getclosestpointto ent pnt))
    (setq ret (CUR:GetCurvature ent pnt))
    (princ "\n离心半径是:")
    (princ ret)
    (and ret (apply 'Ent:Make_Circle ret))
  )
  (princ)
)
;;;=============================================================
;;;测试程序2: 获取椭圆一点处的曲率和离心半径,并比较两种方法
;;;=============================================================
(defun c:tt2 ( / ent obj a b pnt par 1st 2st r1 r2 r3)
  (if (and (setq ent (car (entsel "\n请选取椭圆:")))
	   (setq obj (vlax-ename->vla-object ent))
	   (vlax-property-available-p obj 'MajorRadius)		;这个地方应该加出错处理
  	   (setq a (vla-get-MajorRadius obj))
	   (setq b (vla-get-MinorRadius obj))
      )
    (while (setq pnt (getpoint "\n点取一点:"))
      (setq pnt (trans pnt 1 0))
      (setq pnt (vlax-curve-getclosestpointto ent pnt))
      (setq par (vlax-curve-getparamatpoint ent pnt))
      (setq 1st (vlax-curve-getFirstDeriv ent par))
      (setq 2st (vlax-curve-getSecondDeriv ent par))
      (setq r1  (distance '(0 0) 2st))                          ;这个secondDeriv并不意味着半径
      (setq r2  (cadr (ELL:GetCurvature1 ent a b pnt)))
      (setq r3  (cadr (CUR:getcurvature ent pnt)))
      (princ "\nRadius 1 is:")
      (princ r1)
      (princ "\nRadius 2 is:")
      (princ r2)
      (princ "\nRadius 3 is:")
      (princ r3)
      (UTI:Bench
        10000
        (list
	  (list 'ELL:GetCurvature1 ent a b pnt)
	  (list 'CUR:getcurvature ent pnt)
        )
      )
    )
  )
)
;;;=============================================================
;;;测试程序3: 动态演示曲线的离心半径
;;;=============================================================
(defun c:tt3 (/ CIR ENT LIN PNT PT0 RET 1ST 2ST CEN PAR PT2 VEC)
  (setq ent (car (entsel "\n请选取曲线:")))
  (setq cir (Ent:Make_Circle '(0 0 0) 1))
  (setq lin (ent:make_Line '(0 0 0) '(0 0 0)))
  (setq vec (ent:make_Line '(0 0 0) '(0 0 0)))
  (setq lin (vlax-ename->vla-object lin))
  (setq vec (vlax-ename->vla-object vec))
  (setq cir (vlax-ename->vla-object cir))
  (vla-put-color lin 1)
  (vla-put-color cir 3)
  (vla-put-color vec 6)
  (princ "\n按空格或者回车退出!")
  (while (= (car (setq pt0 (grread 5 0))) 5)
    (setq pnt (trans (cadr pt0) 1 0))
    (setq pnt (vlax-curve-getclosestpointto ent pnt))
    (setq par (vlax-curve-getparamatpoint ent pnt))
    (setq 1st (vlax-curve-getFirstDeriv ent par))
    (setq 2st (vlax-curve-getSecondDeriv ent par))
    (if (setq ret (CUR:GetCurvature ent pnt))
      (progn
	(setq pt2 (mapcar '+ pnt 2st))
	(setq cen (car ret))
        (vlax-put Cir 'Center cen )
        (vlax-put cir 'Radius (cadr ret))
        (vlax-put lin 'StartPoint cen)
        (vlax-put lin 'EndPoint pnt)
	(vlax-put vec 'StartPoint pt2)
        (vlax-put vec 'EndPoint pnt)
      )
    )
  )
  (vla-erase cir)
  (vla-erase lin)
  (vla-erase vec)
  (princ)
)
;;;=============================================================
;;;测试程序4: 由曲线的离心半径描绘其轨迹
;;;=============================================================
(defun c:tt4 (/ ent lst par px1 px2 pxn pnt Inf cen)
  (if (setq ent (car (entsel "\n请选取曲线:")))
    (progn
      (setq px1 (vlax-curve-getStartParam ent))
      (setq px2 (vlax-curve-getendparam ent))
      (setq pxn (* (- px2 px1) 0.02))
      (setq par px1)
      (while (<= par px2)
        (setq pnt (vlax-curve-getpointatparam ent par))
        (setq Inf (CUR:GetCurvature ent pnt))
        (if (setq cen (car Inf))
	  (setq lst (cons cen lst))
	)
        (setq par (+ par pxn))
      )
      (setq lst (reverse lst))
      (Ent:Make_LWPoly lst 1)
    )
  )
)

相关联代码:

;;;************************************************************;
;;;实体创建部分
;;;************************************************************;
;;;-------------------------------------------------------------
;;;创建一个点
;;;输入: 一个三维或者二维的点
;;;输出: 点实体的图元名
;;;-------------------------------------------------------------
(defun Ent:Make_Point (p)
  (entmakex (list '(0 . "POINT") (cons 10 p)))
)
;;;-------------------------------------------------------------
;;;创建一个带颜色的点(此函数为测试或者其他用途)
;;;输入: 一个三维或者二维的点表和一个颜色号
;;;输出: 点实体的图元名
;;;-------------------------------------------------------------
(defun Ent:MakePoint-1 (p c)
  (entmakex (list '(0 . "POINT") (cons 10 p) (cons 62 c)))
)
;;;-------------------------------------------------------------
;;;创建一条直线段
;;;输入: 两个三维或者二维的点
;;;输出: 线段实体的图元名
;;;-------------------------------------------------------------
(defun Ent:Make_Line (p q)
  (entmakeX (list '(0 . "LINE") (cons 10 p) (cons 11 q)))
)
;;;-------------------------------------------------------------
;;;创建一个由三条直线组成的三角形
;;;输入: 三个三维或者二维的点
;;;输出: 由三条直线组成的三角形
;;;-------------------------------------------------------------
(defun Ent:Make_Triangle (p1 p2 p3)
  (mapcar 'Ent:Make_Line (list p1 p2 p3) (list p2 p3 p1))
)
;;;-------------------------------------------------------------
;;;创建一个三维多段线
;;;输入: 三维的点集
;;;输出: 三维多段线实体
;;;-------------------------------------------------------------
(defun Ent:Make_Poly (pts Closed / e)
  (if Closed
    (setq e (Entmake (list '(0 . "POLYLINE") '(70 . 9))))
    (setq e (Entmake (list '(0 . "POLYLINE") '(70 . 8))))
  )
  (foreach p pts
    (entmake (list '(0 . "VERTEX") '(70 . 32) (cons 10 p)))
  )
  (entmake '((0 . "SEQEND")))
  (entlast)
)
;;;-------------------------------------------------------------
;;;创建轻多段线
;;;输入: 二维的点集
;;;输出: 轻多段线实体名
;;;-------------------------------------------------------------
(defun Ent:Make_LWPoly (pts closed /)
  (entmakeX
    (append
      '((0 . "LWPOLYLINE")
        (100 . "AcDbEntity")
        (100 . "AcDbPolyline")
       )
      (list (cons 90 (length pts)))                      	;顶点个数
      (mapcar (function (lambda (x) (cons 10 x))) pts)  	;多段线顶点
      (list (cons 70  closed))                  	        ;闭合的
    )
  )
)
;;;-------------------------------------------------------------
;;; 创建圆实体
;;; 输入: 圆心C和半径R
;;; 输出: 圆的实体名
;;;-------------------------------------------------------------
(defun Ent:Make_Circle (C R)
  (entmakex (list '(0 . "CIRCLE") (cons 10 C) (cons 40 R)))
)

程序1效果演示:

程序2效果演示:

LISP的表达式求值

表达式一般来说有三种:前缀表达式、中缀表达式、后缀表达式,其中后缀表达式又叫做逆波兰表达式。中缀表达式是最符合人们思维方式的一种表达式,顾名思义,就是操作符在操作数的中间。而前缀表达式和后缀表达式中操作符分别在操作数的前面和操作数的后面。在写表达式,我们一般用中缀表达式,譬如 1+2*3-4/5。并且按照操作符的优先级进行计算。
然而LISP语言是一种前缀表达式,为了把表达式转为LISP函数或者求值,需要进行翻译,添加大量的括号和修改函数的顺序。
这个程序的目的就是使得这一工作变简单。
当然,CAD里面本身也有几种种方式能完成这个,但它们的优缺点容我后面讨论。
程序借鉴了飞诗的一些代码,在此深表感谢。
程序的核心代码如下:

;;;=============================================================
;;; 函数目的: 字符表达式转为函数,主要用于多次调用时提升速度
;;; 输入: expr--字符表达式,sFunc--函数名,sArg--参数符号列表
;;; 输出: 定义函数,并返回其名
;;; 例子: (CAL:Expr2Func "sin(x)+20*y" 'test '(x y))
;;; 结果: 定义了一个名为test的函数,参数符号为x y
;;; 注意: 除法区分整数和浮点数,譬如"2/3"和"2/3.0"结果不同;
;;;       可用自定义函数,前提是首先要加载;
;;;       可用科学计算法,但应满足LISP中的语法。建议用括号;
;;;       表达式应满足语法要求,小数和乘号不能按省略写法。
;;;=============================================================
(defun CAL:Expr2Func (expr sFunc sArgs / lst)
  (setq lst (CAL:Separate expr))				;先按照括号空格和运算符分离字符
  (setq lst (CAL:Operators lst '((^ . expt)) ()))	        ;乘方(幂)最优先
  (setq lst (CAL:Operators lst '((* . *) (/ . /) (% . rem)) ()));其次乘除和求模运算
  (setq lst (CAL:Operators lst '((+ . +) (- . -)) ()))		;最后处理加减法运算
  (defun-q-list-set sFunc (cons sArgs lst))			;表达成函数
  sFunc
)
;;;=============================================================
;;; 函数目的: 字符表达式求值
;;; 输入: expr--字符表达式
;;; 输出: 计算表达式的结果
;;; 例子: (CAL:Expr2Value "sin(1)+20*2")
;;; 结果: 40.8415
;;;=============================================================
(defun CAL:Expr2Value (expr / lst)
  (setq lst (CAL:Separate expr))				;先按照括号空格和运算符分离字符
  (setq lst (CAL:Operators lst '((^ . expt)) ()))	        ;乘方(幂)最优先
  (setq lst (CAL:Operators lst '((* . *) (/ . /) (% . rem)) ()));其次乘除和求模运算
  (setq lst (CAL:Operators lst '((+ . +) (- . -)) ()))		;最后处理加减法运算
  (eval (car lst))						;求值
)
;;;=============================================================
;;; 函数目的: 先分离出函数和+-*/%^运算符,其余均视作变量或数值,
;;; 并简单检查括号匹配。
;;; 输入: expr--字符表达式
;;; 输出: 函数(包括运算符)和变量及数值的列表
;;;=============================================================
(defun CAL:Separate (expr / CHAR FUNS LASTCHAR LST Temp LBRACKET RBRACKET next)
  (setq expr (vl-string-translate "{[]}\t\n," "(())   " expr))  ;替换{[]}\t\n,字符
  (setq expr (strcase expr t))					;全部转为小写
  (setq funs '("+" "-" "*" "/" "^" "%" ))		        ;按照基本运算符分割字符
  (setq Temp "")
  (setq lst "(")
  (setq Lbracket 0)						;左括号计数器
  (setq Rbracket 0)						;右括号计数器
  (while (/= expr "")
    (setq char (substr expr 1 1))                               ;字符串的第一个字符
    (setq next (substr expr 2 1))				;字符串的第二个字符
    (if	(or (= char "(")
	    (= char ")")					;括号一定是分隔符
	    (and (= char " ") (/= next "(") (/= next " "))      ;如果不是连续的空格符
	    (and (member char funs)				;根据运算符进行分割
	         (not (CAL:isScientific temp lastchar char))    ;忽略科学计数法
	    )
	)
      (progn
	(if (CAL:IsFunction (Read temp))			;如果为普通函数
	  (setq	lst	 (strcat lst "(" Temp " " )		;则把括号移至函数符号前
		Lbracket (1+ Lbracket)				;左括号计数器加1
	  )
	  (progn
	    (and (= char "(") (setq Lbracket (1+ Lbracket)))    ;左括号计数器加1
	    (and (= char ")") (setq Rbracket (1+ Rbracket)))	;右括号计数器加1
	    (setq lst (strcat lst Temp " " char " "))
	  )
	)
	(setq Temp "")                                          ;如果是函数或者括号空格之类,则在此处重新开始
      )
      (setq Temp (strcat Temp char))                            ;否则连取这个字符
    )
    (setq expr (substr expr 2))					;字符串剩下的字符
    (setq lastchar char)
  )
  (if (/= Lbracket Rbracket)					;如果括号不平衡
    (alert "括号不匹配(Mismatched Brackets)!")			;警告信息
    (read (strcat lst Temp ")"))				;否则转为表
  )
)
;;;=============================================================
;;; 函数目的: 分析+-*/%^运算符,并组合到表中
;;; 输入: lst-已分割的表,funs-待分析的运算符,Recursive-是否递归
;;; 输出: 函数(包括运算符)和变量及数值的列表
;;;=============================================================
(defun CAL:Operators (lst funs Recursive / fun L n)
  (foreach a lst
    (if	(listp a)
      (setq a (CAL:Operators a funs T))				;如果元素为表,则递归进去
    )
    (if (= (type a) 'INT)
      (setq a (float a))
    )
    (if	(setq fun (cdr (assoc (car L) funs)))                   ;前一个符号为+-*/%^运算符
      (if (or (null (setq n (cadr L)))                          ;前前一个符号为空
	      (and (VL-SYMBOLP n) (CAL:IsFunction n))           ;或者是函数符号
	  )
	(setq L (cons (list fun a) (cdr L)))                    ;无须交换位置
	(setq L (cons (list fun n a) (cddr L)))	                ;交换运算符和操作数位置
      )
      (setq L (cons a L))                                       ;其他的不做改变
    )
  )
  (setq n (car L))
  (if (and Recursive (not (cadr L)) (or (listp n) (numberp n))) ;如果是递归的,而且只有一个元素,且这个元素为表或者数字
    n								;那么就只取这个元素,以防止多余括号出现
    (reverse L)							;cons运算后的反转表列
  )
)
;;;=============================================================
;;; 函数目的: 判断一个符号是不是普通函数(内部函数或自定义函数)
;;;=============================================================
(defun CAL:IsFunction (n / s)
  (setq s (type (eval n)))
  (and (or (= s 'SUBR) (= s 'USUBR)) (not (member n '(+ - * /))))
)
;;;=============================================================
;;; 函数目的: 检测一个字符串是否是科学计数法(是否有更好方法?)
;;;=============================================================
(defun CAL:isScientific (temp lastchar char)
  (and (= lastchar "e") (numberp (read (strcat temp char "0"))))
)
;;;=============================================================
;;; 函数目的: 检查函数表达式转函数的结果
;;; 输入: lst,用cal:expr2func求得的表
;;; 输出: 如果表达式里有非参数且未赋值的变量符号则返回nil, 否则T
;;; 例子: (CAL:CheckFunc (CAL:Expr2func "sin(a)+20*2" 'fx '(x)))
;;; 结果: nil
;;;=============================================================
(defun CAL:CheckFunc (lst / isOK CAL:TempSym Args)
  (setq IsOK T)
  (setq Args (car lst))
  (while (setq lst (cdr lst))
    (setq CAL:TempSym (car lst))                                ;对表中的每个元素
    (if	(listp CAL:TempSym)					;如果这个元素为表
      (if CAL:TempSym
	(setq IsOk (CAL:CheckFunc (cons Args CAL:TempSym)))	;且不为空则递归进去
	(setq IsOk nil)                                         ;否则检测结果为假
      )
      (if (and (vl-symbolp CAL:TempSym)                         ;如果是一个符号
	       (not (member CAL:TempSym Args))			;且不为参数表中的符号
	       (not (vl-symbol-value CAL:TempSym))              ;且未赋值
	  )
	(setq IsOk nil)						;则检测结果为假
      )
    )
    (if	(null IsOK)
      (setq lst nil)
    )
  )
  IsOK
)

;;;=============================================================
;;;以下函数为自定义的一些简单的数学函数
;;;=============================================================
(defun r2d (x) (* 57.2957795130823208768 x))                    ;弧度转度
(defun d2r (x) (* 0.01745329251994329577 x))                    ;度转弧度
(defun int (x) (atoi (rtos x 2 0)))                             ;四舍五入取整函数
(defun ceil (x) (1+ (fix x)))                                   ;天花板函数
(defun ln (x) (log x))            				;以e为底的对数函数
(defun log10 (x) (* (log x) 0.43429448190325182765))            ;以10为底的对数函数
(defun exp10 (x) (expt 10 x))					;以10为底的指数函数
(defun pow (x y) (expt x y))                                    ;指数函数
(defun tan (x) (/ (sin x) (cos x)))				;正切函数
(defun cot (x) (/ (cos x) (sin x)))				;余切函数
(defun sec (x) (/ 1 (cos x)))                                   ;正割函数
(defun csc (x) (/ 1 (sin x)))					;余割函数
(defun asin (x) (atan x (sqrt (- 1 (* x x)))))                  ;反正弦函数
(defun acos (x) (atan (sqrt (- 1 (* x x))) x))			;反余弦函数
(defun sinh (x) (* 0.5 (- (exp x) (exp (- x)))))		;双曲正弦函数
(defun cosh (x) (* 0.5 (+ (exp x) (exp (- x)))))		;双曲余弦函数
(defun tanh (x) (- 1 (/ 2 (1+ (exp (+ x x)))))) 		;双曲正切函数
(defun coth (x) (/ 1 (tanh x)))					;双曲余切函数
(defun sech (x) (/ 1 (cosh x)))					;双曲正割函数
(defun csch (x) (/ 1 (sinh x)))					;双曲余割函数
(defun asinh (x) (log (+ x (sqrt (1+ (* x x))))))		;反双曲正弦函数=log(x+sqrt(x*x+1))
(defun acosh (x) (log (+ x (sqrt (1- (* x x))))))       	;反双曲余弦函数=log(x+sqrt(x*x-1))
(defun atanh (x) (log (sqrt (/ (+ 1 x)(- 1 x)))))		;反双曲正切函数=log(sqrt((1+x)/(1-x)))
(defun revSign (x) (- x))					;反号函数
(defun reciprocal (x) (/ 1.0 x))				;倒数
(defun sqr (x) (* x x))						;平方函数
(defun cube (x) (* x x x))					;立方函数
(defun cuberoot	(x)						;立方根函数
  (if (minusp x)
    (- (expt (- x) 0.333333333333333333333))
    (expt x 0.333333333333333333333)
  )
)
(defun round (x / y)						;四舍五入函数
  (setq y (fix x))
  (if (< (abs (- x y)) 0.5)
    y
    (if (< x 0)
      (1- y)
      (1+ y)
    )
  )
)

以下是一些测试:

;;; 例子:
;;; (CAL:Separate "(sin(-x)-cos(-x+(1+8*(2/7))+2^4-5))*0.5-0.5e-20+20*cos(x)+20")
;;; 结果: ((SIN - X) - (COS - X + (1 + 8 * (2 / 7)) + 2 ^ 4 - 5))
;;; (CAL:Expr2Func "(sin(+x)-cos(-x+(1+8*(2/7))+(2^4)-5))*0.5-0.5e-20+20*cos(x)+20" 'test '(x))
;;; 结果: 定义了一个名为test的函数,参数符号为x
;;; (CAL:Expr2Value "(sin(+0.5)-cos(-pi+(1+8*(2/7))+(2^4)-5))*0.5-0.5e-20+20*cos(pi/2)+20")
;;; 结果: 20.6616
;;; 以下是关于这个程序的其他方法:
;;; 方法一:用cal函数计算
;;; 如:(cal "1+4+5*2+(5+5)/2+((6+6)/2+(5+5)/2)")
;;; 优点:CAD内置函数。
;;; 缺点:这个函数要求先要加载cal函数.并且三角函数会自动把变量或者数值理解为角度。
;;; 方法二:wcs脚本语言法,无痕提出的一种方法
;;; (setq wcs (vla-GetInterfaceObject (vlax-get-acad-object) "ScriptControl"))
;;; (vlax-put-property wcs "language" "vbs")
;;; (vla-eval wcs "1+4+5*2+(5+5)/2+((6+6)/2+(5+5)/2)")  ;返回 ->31.0
;;; 优点:能按照vb的语法直接计算。
;;; 缺点:难以定义表达式为函数,不能利用自定义函数,在64位CAD上此法行不通,因为不能创建脚本对象。
;;; 下面例子为在CAD中绘制函数图像
(defun c:test1(/ expr a b d x y e pts)
  (setq expr (getstring "\n请输入表达式:"))
  (initget 1)
  (setq a (getreal "\n上届:"))
  (initget 1)
  (setq b (getreal "\n下届:"))
  (if (CAL:EXPR2FUNC  expr 'test '(x))
    (progn
      (setq d (/ (- b a) 1000.0))
      (setq x a)
      (setq pts nil)
      (repeat 1000
	(setq x (+ x d))
	(setq y (test x))
	(setq pts (cons (list x y 0) pts))
      )
      (setq pts (reverse pts))
      (setq e (Entmake (list '(0 . "POLYLINE") '(70 . 8))))
      (foreach p pts
        (entmake (list '(0 . "VERTEX") '(70 . 32) (cons 10 p)))
      )
      (entmake '((0 . "SEQEND")))
      (entlast)
    )
  )
)
;;; 在CAD中绘制参数曲线
;;; x=a*(2*cos(t)-cos(2*t))
;;; y=a*(2*sin(t)-sin(2*t))
(defun c:test2 (/ expr1 expr2 a b d k x y pts e)
  (setq expr1 "3*(2*cos(k)-cos(2*k))")
  (setq expr2 "3*(2*sin(k)-sin(2*k))")
  (setq a 0)
  (setq b (+ pi pi))
  (CAL:EXPR2FUNC expr1 'fx '(k))
  (CAL:EXPR2FUNC expr2 'fy '(k))
  (setq d (/ (- b a) 360))
  (setq k a)
  (setq pts nil)
  (repeat 360
    (setq k (+ k d))
    (setq x (fx k))
    (setq y (fy k))
    (setq pts (cons (list x y 0) pts))
  )
  (setq pts (reverse pts))
  (setq e (Entmake (list '(0 . "POLYLINE") '(70 . 9))))
  (foreach p pts
    (entmake (list '(0 . "VERTEX") '(70 . 32) (cons 10 p)))
  )
  (entmake '((0 . "SEQEND")))
  (entlast)
)
;;; 定义为函数后,明显速度快多了
(defun c:test3(/ str1 str2 x)
  (setq str1 "(sin(+x)-cos(-x+(1+8*(2/7.0))+(2^4)-5))*0.5-0.5e-20+20*cos(x)+20")
  (setq str2 "(sin(r2d(x))-cos(r2d(-x+(1+8*(2/7.0))+(2^4)-5)))*0.5-0.5e-20+20*cos(r2d(x))+20")
  (CAL:Expr2Func str1 'f1 '(x))
  (setq x 12)
  (uti:bench 1000
    (list
      (cons 'f1 '(12))
      (cons 'CAL:Expr2Value (list str1))
      (cons 'cal (list str2))
    )
  )
)