王依依

王依依的博客

她的个人主页  她的博客

对于lisp党来说, 一切控制结构都是纸老虎

王依依  2009年10月19日 星期一 10:44 | 3692次浏览 | 96条评论

数值分析, 非线性方程求解, 做作业用的.

;; -*- coding:utf-8; syntax:common-lisp -*-

(defun bisection-recur (f a b &optional (e 1.0L-6) (recur-nth 0))
  "Page 79. bisection"
  (let ((alpha (/ (+ a b) 2.0L0)))            ; (1) use long-float
    (format t "k=~D a=~F b=~F f(alpha)=~F |b-a|=~F~%"
            recur-nth a b (funcall f alpha) (abs (- b a)))
    (cond ((or (< (abs (- b alpha)) e) 
               (< (abs (funcall f alpha)) e))
            alpha)                      ; (2)
          ((< (* (funcall f alpha) 
                 (funcall f a)) 
              0)
           (bisection f a alpha e (1+ recur-nth))) ;(3)
          (t (bisection f alpha b e (1+ recur-nth)))))) ; (4)
          

(defun bisection (f a b &optional (e 1.0L-6)
                   &key (max-recur 100))
  "Page 79. bisection"
  (loop 
     :for k from 0 to max-recur
     :for alpha = (/ (+ a b) 2.0L0)     ; (1)
     :do (format t "k=~D a=~F b=~F f(alpha)=~F |b-a|=~F~%"
                 k a b (funcall f alpha) (abs (- b a)))
     :thereis (when (or (< (abs (- b alpha)) e) 
                        (< (abs (funcall f alpha)) e))
                alpha)                  ; (2)
     :when (< (* (funcall f alpha) (funcall f a)) 0) :do (setq b alpha) ; (3)
     :else :do (setq a alpha)))         ; (4)

;; TEST: (simple-iteration-aitken #'testphi 5)
(defun simple-iteration-aitken (phi x0 &optional (e 1.0L-6)
                                    &key (max-recur 100))
  "Page 87. Aitken speedup."
  (let ((x (coerce x0 'long-float))
        (alpha nil))
    (format t "x_~D=~F~%" 0 x)
    (loop
       :for k from 0 to max-recur
       :for y = (funcall phi x)
       :for z = (funcall phi y)
       :when (= 0 (+ (- z (* 2 y)) x)) :do (return x)
       :do (setq alpha (- x
                          (/ (expt (- y x) 2)
                             (+ (- z (* 2 y)) x))))
       :do (format t "x_~D=~F |x_~D-x_~D|=~F~%"
                   (1+ k) alpha (1+ k) k (abs (- alpha x)))
       :when (< (abs (- alpha x)) e) :do (return alpha)
       :do (setq x alpha))))
     

;; TEST: (newton-iteration #'testfun #'testfun^ 5)
(defun newton-iteration (f f^ x0 &optional (e 1.0L-6)
                           &key (max-recur 100))
  "Page 92. Newton iter."
  (let ((x (coerce x0 'long-float))
        (alpha nil))
    (loop 
       :for k from 0 to max-recur
       :for fx = (funcall f x)
       :for f^x = (funcall f^ x)        ; (2)
       :do (format t "k=~D x_~D=~F f(x_~D)=~F "
                   k k x k fx)
       :when (= f^x 0) :do (return nil) ; (3)
       :do (setq alpha (- x (/ fx f^x))) ; (4)
           (format t "|x_~D-x_~D|=~F ~%" 
                   (1+ k) k (abs (- alpha x)))
       :thereis (when (< (abs (- alpha x)) e) ; (5)
                  alpha)
       :do (setq x alpha))))            ; (7)


;; TEST: (newton-iteration-simplify #'sin 0.7 5)
(defun newton-iteration-simplify (f M x0 &optional (e 1.0L-6)
                                    &key (max-recur 100))
  "Page 92. Newton iter simplify."
  (let ((x (coerce x0 'long-float))
        (alpha nil))
    (loop 
       :for k from 0 to max-recur
       :for fx = (funcall f x)
       :do (format t "k=~D x_~D=~F f(x_~D)=~F "
                   k k x k fx)
           (setq alpha (- x (/ fx M)))
           (format t "|x_~D-x_~D|=~F ~%" 
                   (1+ k) k (abs (- alpha x)))
       :thereis (when (< (abs (- alpha x)) e)
                  alpha)
       :do (setq x alpha))))

;; TEST: (secant-method #'testfun 5 4)
(defun secant-method (f x0 x1 &optional (e 1.0L-6)
                        &key (max-recur 100)x)
  "Page 92. Newton iter simplify."
  (let ((xk_1 (coerce x0 'long-float))
        (xk (coerce x1 'long-float))
        (alpha nil))
    (loop 
       :for k from 1 to max-recur
       :for fx = (funcall f xk)
       :for f^x = (/ (- fx (funcall f xk_1))
                     (- xk xk_1))
       :do (format t "k=~D x_~D=~F f(x_~D)=~F "
                   k k xk k fx)
       :when (= f^x 0) :do (return nil)
       :do (setq alpha (- xk (/ fx f^x)))
           (format t "|x_~D-x_~D|=~F ~%" 
                   (1+ k) k (abs (- alpha xk)))
       :thereis (when (< (abs (- alpha xk)) e)
                  alpha)
       :do (setq xk_1 xk)
           (setq xk alpha))))




;; test case
        
(defun testfun (x)
  (- (* x (exp x)) 1))

(defun testfun^ (x)
  (+ (exp x)
     (* x (exp x))))

;; use alpha=pi/2 to test
(defun testphi (x)
  (+ 1.6 (* 0.99 (cos x))))

评论

我的评论:

发表评论

请 登录 后发表评论。还没有在Zeuux哲思注册吗?现在 注册 !
李书冰

回复 李书冰  2010年04月07日 星期三 15:17

围观片刻…掩面叹息…深受打击的80后……

0条回复

Shawn

回复 Shawn  2009年10月26日 星期一 11:49

久违的括号...这是common lisp还是?lisp语言只用过scheme~建议下次帖code前把问题的数学模型也贴上来不然简直看不懂

0条回复

劳永超

回复 劳永超  2009年10月23日 星期五 00:00

很牛~ 我们以前数值分析还不用写代码的..

0条回复

胡锦涛

回复 胡锦涛  2009年10月22日 星期四 16:38

又看见可爱的括号哈哈~
只研究EL~

4条回复

  • 王依依

    回复 王依依  2009年10月22日 星期四 16:40

    同是 emacs 党~哈哈.

    3条回复

      • 胡锦涛

        回复 胡锦涛  2009年10月22日 星期四 16:46

        还是emacs buffer方便啊,私活 小游戏 email 什么都不耽误……

        2条回复

          • 王依依

            回复 王依依  2009年10月22日 星期四 16:48

            有没有能把 mplayer 嵌进去的 mode?

            1条回复

              • 胡锦涛

                回复 胡锦涛  2009年10月22日 星期四 16:57

                可能还是有的,不过估计会为此放弃太多.
                为了树立办公室榜样形象 我还是放弃变为耳机党了……
                本来emacs启动就跟爬一样,再加个mplayer进去.就能一大清早来开个emacs然后倒茶 看报纸 查询违章记录再上个厕所,回来就开了.
                如buffer切换牛速就是不能忍的~

                0条回复

赵斌

回复 赵斌  2009年10月21日 星期三 23:40

这只猫不一般....
一般的猫不喜欢,这只猫得拜一

0条回复

李波

回复 李波  2009年10月21日 星期三 10:30

很好很强大

0条回复

龙黎

回复 龙黎  2009年10月21日 星期三 08:39

好强大的猫猫,Lisp 党也可叫“括号党”吧 :)

1条回复

徐继哲

回复 徐继哲  2009年10月21日 星期三 00:13

这篇博客这么火爆啊。

2条回复

alexpeng

回复 alexpeng  2009年10月20日 星期二 22:21

支持猫猫.
顺便问一句,你的老师用Lisp吗?

1条回复

猫哥

回复 猫哥  2009年10月20日 星期二 22:05

呃~你有没想过写一高阶过程自动识别传入的方程式?而不是每次为一个方程设计一个过程
另外2分为什么要写两个实现?还是就是为了实现而已

6条回复

  • 王依依

    回复 王依依  2009年10月22日 星期四 11:23

    识别方程式涉及到 parse 表达式解析, 这方面倒不如用 haskell 的 GADT方便

    2条回复

      • 猫哥

        回复 猫哥  2009年10月24日 星期六 12:45

        呵呵~要是只图方便,那都用matlab好了~用基础语言就是为了写着玩嘛

        1条回复

  • 王依依

    回复 王依依  2009年10月22日 星期四 11:20

    自动传入的话参数个数不一样, 另外也没多少重复方法
    最多可以把牛顿迭代法及其简化方法做一个抽象, 意义不大

    二分写着玩, 第一个是尾递归, 第二个用 loop

    2条回复

      • 猫哥

        回复 猫哥  2009年10月24日 星期六 12:48

        好像你没理解我的意思,不过也无所谓了~参数个数不一样跟我说的那个没太大关联,不过这都是题外话了,祝你玩得愉快哈

        1条回复

          • 王依依

            回复 王依依  2009年10月24日 星期六 15:34

            知道, 如果几和函数里有很类似的代码段, 那么很可能有个宏抽象可以处理这些.
            这个我想过, 不太好做.

            祝你也玩得愉快.

            0条回复

Shadow

回复 Shadow  2009年10月20日 星期二 20:01

专程来膜拜的or2……

0条回复

张中原

回复 张中原  2009年10月20日 星期二 19:46


实际中用到这个语言的几率很低很低
种语言本身都有自己的思维方式
工程和学术

0条回复

電波系山寨文化科学家

回复 電波系山寨文化科学家  2009年10月20日 星期二 19:36

曾经学会过lisp

0条回复

华雪亮

回复 华雪亮  2009年10月20日 星期二 19:21

用Scheme来弄读起来会简单点..... :)

4条回复

  • 王依依

    回复 王依依  2009年10月22日 星期四 11:20

    scheme 表达能力弱很多

    3条回复

      • 华雪亮

        回复 华雪亮  2009年10月22日 星期四 11:50

        呵呵...Plt Scheme ...应该表达你这个不成问题吧....没怎么用过Common Lisp..elisp对我来说足够了...

        2条回复

          • 王依依

            回复 王依依  2009年10月22日 星期四 16:47

            嗯, plt 很强大的说.
            那我应该说.....
            是代码信息量吧....

            1条回复

              • 华雪亮

                回复 华雪亮  2009年10月22日 星期四 17:02

                不过能把lisp写成那样,真不简单那...值得学习学习.....

                0条回复

端瑞

回复 端瑞  2009年10月20日 星期二 15:09

这个不顶不行啊!加油!

0条回复

于昕

回复 于昕  2009年10月20日 星期二 14:13

标题真强悍啊.

0条回复

李登

回复 李登  2009年10月20日 星期二 14:11

强~

0条回复

Mr.Gui

回复 Mr.Gui  2009年10月20日 星期二 13:31

MM太厉害了!! lisp都会。
我现在都27了, 才刚刚接触编程, 而且还是大家一致公认简单的python, 可就是python我学得还是磕磕绊绊。 只能佩服你们了。 各位都是怎么炼成的, 有什么比较科学的学习方法吗?

1条回复

  • 王依依

    回复 王依依  2009年10月22日 星期四 11:21

    python 不错啊~ 很好学的
    简明看完然后看深入点的教材

    0条回复

李建

回复 李建  2009年10月20日 星期二 12:58

要想倒腾计算机更有趣点(或者说有水平点),Lisp 真是跳不过的一道关,最近又拿起 Emacs Lisp Info 在看。LZ 的程序,我现在水平还看不懂,唯有顶贴了。

0条回复

杨昆

回复 杨昆  2009年10月20日 星期二 12:11

mm是哪个学校的阿,居然学lisp.听说我们学校的计算机入门语言改为额schema了,难不成是校友

1条回复

暂时没有评论

Zeuux © 2024

京ICP备05028076号