maximaの計算時間を計測

maximaで

showtime:true;

とすれば、以後計算するたびに計算時間が表示される。

meval

defmfun で定義されたlispの関数はmevalで評価できる。

(defmfun $test (&rest keys)
  (let ()
    (setf keys (push '$factorout keys))
    (meval keys)
    ))

自作関数のメモ(プロトタイプの時の)

式に指定した関数が使われているか調べる。式に積分、伊藤積分が含まれる場合、その要素の中も調べる。

depends_func (expr func)
> true (if expr depends func)
> false (otherwise)
ex.
(depends_func expr 'mplus)

file:

file:

file:

lisp内で記述するmaximaの数学記号

自然対数

'$%E

maximaで使用する数式の関数f(x):=x+1などをlisp内で定義する

(setf (symbol-plist $f) '(MPROPS (NIL MEXPR ((LAMBDA) ((MLIST) $x) (MPLUS $x 1)))))

maximaで
f(x):=x+1;
としたことと同義。

下記のようにすれば、maximaの表示上は定義してように見えるが、その関数を呼び出すことはできない。

(list `(mdefine simp) func2 (mul2 ini (pow '%e ($intstyle e2 y))))

メモ

(setq func1 (concat x n))
(setq func2 (list `(,func1 simp) y))
(setf result1 (mul2 ini (pow '%e ($intstyle e2 y))) )
(setf result2 `(MPROPS (NIL MEXPR ((LAMBDA) ((MLIST) ,y) ,result1))) )
(setf (symbol-plist func1) result2)

maxima <=> lisp

f(s):=s^2
((MDEFINE SIMP) (($F) $S) ((MEXPT) $S 2)) 

fubiniの定理

確率積分と時間微分の交換はできた。

今後、このページに課題を書く。

assume

maximaを使っていて、
「positive or negative?」
と返されるが、assume(a > 0) とかとしておけばオーケー

Fubiniの定理

Fubiniの定理の実装

目標:時間積分と確率積分の順序を交換。

task:
1. 式からintegrate or ito_integrate を探す。

2. integrandからito_integrate または integrate を探す。
2-1. integrate(ito_integrate())となっているのを発見 => 3
2-2. otherwise => 2

3.順序交換

式の展開 ratsimp ratf ratdisrep

file: rat3e.lisp
ん~ よくわからん。とりあえず、ratsimpをすればいいんだな。

(defun ratsimp (x varlist genvar) ($ratdisrep (ratf x)))

(defmfun ratf (l)
  (prog (u *withinratf*)
     (setq *withinratf* t)
     (when (eq '%% (catch 'ratf (newvar l)))
       (setq *withinratf* nil)
       (return (srf l)))
     (setq u (catch 'ratf (ratrep* l))) ; for truncation routines
     (return (or u (prog2 (setq *withinratf* nil) (srf l))))))

(defmfun ratdisrep (e) (simplifya ($ratdisrep e) nil))

depends

file: comm.lisp
eがxに依存しているか判別する。

(defun depends (e x)
  (cond ((alike1 e x) t)
 ((mnump e) nil)
 ((atom e) (mget e 'depends))
 (t (or (depends (caar e) x) (dependsl (cdr e) x)))))

伊藤の公式

伊藤の公式の実装。
ブラウン運動が異なるときはちゃんと二次変分がゼロになってます。

問題点
1. 初期値を考慮するのを忘れてる。FINISHED!!
2. 時間微分を忘れている。FINISHED!!


課題
1. 被積分関数の整理。係数を前に出す、項を分けるなど。FINISHED 10/11/09
2. \int (\int f dw \int g dw) dw 中も伊藤の公式で整理する。FINISHED!!


発展
1. Fubiniの定理
2. 確率微分方程式

diffint, diffint1

file: comm2.lisp
diff(integrate(f(s,x),s,t,x),x)を計算する過程で呼ばれる。
不定積分=>
定積分=> 積分と微分の交換、積分範囲で微分

(defmfun diffint (e x) ;; example ((%INTEGRATE . #1=(SIMP)) (($F . #1#) $S) $S 0 $T) $X
  (let (a)
    (cond 
      ;;indefinite integral
      ((null (cdddr e)) 
           (cond ((alike1 x (caddr e)) (cadr e))
             ((and (not (atom (caddr e))) (atom x) (not (free (caddr e) x)))
              (mul2 (cadr e) (sdiff (caddr e) x)))
             ((or ($constantp (setq a (sdiff (cadr e) x)))
                  (and (atom (caddr e)) (free a (caddr e))))
              (mul2 a (caddr e)))
             (t (simplifya (list '(%integrate) a (caddr e)) t))))
      ;;definite integral
      ;;integrating variable = differential variable 
      ((alike1 x (caddr e)) (addn (diffint1 (cdr e) x x) t))
      ;;integrating variable =/ differential variable 
      (t (print (simplifya (list 'b) t)))
      (t (addn (cons 
                     ;; differentiate integrand with respect to x
                     (if (equal (setq a (sdiff (cadr e) x)) 0)
                         0
                         ;; integrate a 
                         (simplifya (list '(%integrate) a (caddr e) (cadddr e) (car (cddddr e))) t))
                     ;; arguments are (integrand, integrating variable, integral range), differential variable, integrating variable
                     (diffint1 (cdr e) x (caddr e))) t)))))


(defun diffint1 (e x y)
  (let ((u (sdiff (cadddr e) x)) (v (sdiff (caddr e) x)))
    ;;int_(caddr e)^(cadddr e) (car e)(y) dy
    ;;(mul2 'a 'b) => '((%mtimes) a b)
    (list (if (pzerop u) 0 (mul2 u (maxima-substitute (cadddr e) y (car e))))
   (if (pzerop v) 0 (mul3 v (maxima-substitute (caddr e) y (car e)) -1)))))

addn

file: opers.lisp
渡されたlistの中身をすべて足す。

(defmfun addn (terms simp-flag)
  (cond ((null terms) 0)
         (t (simplifya `((mplus) . ,terms) simp-flag))))

nonvarcheck

file: comm.lisp

(defmfun nonvarcheck (e fn)
  (if (or (mnump e)
   (maxima-integerp e)
   (and (not (atom e)) (not (eq (caar e) 'mqapply)) (mopp1 (caar e))))
      (merror (intl:gettext "~:M: second argument must be a variable; found ~M") fn e)))

deriv

file:comm.lisp

(defun deriv (e) ;; e <- (((MEXPT SIMP) $X 2) $X 1)
  (prog (exp z count) 
           (cond ((null e) (wna-err '$diff))
    ((null (cdr e)) (return (stotaldiff (car e))))
    ((null (cddr e)) (nconc e '(1))))
     (setq exp (car e) z (setq e (copy-list e)))
     loop (if (or (null derivlist) (member (cadr z) derivlist :test #'equal)) (go doit))
     ; DERIVLIST is set by $EV
     (setq z (cdr z))
     loop2(cond ((cdr z) (go loop))
  ((null (cdr e)) (return exp))
  (t (go noun)))
     doit (cond ((nonvarcheck (cadr z) '$diff))
  ((null (cddr z)) (wna-err '$diff)) ;;(cddr z) <- (1)
  ((not (eq (ml-typep (caddr z)) 'fixnum)) (go noun)) ;;(caddr z) <- 1
  ((minusp (setq count (caddr z))) ;; count <- 1
   (merror (intl:gettext "diff: order of derivative must be a nonnegative integer; found ~M") count)))
     ;;z <-((MEXPT SIMP) $X 2)
     loop1 (cond ((zerop count) (rplacd z (cdddr z)) (go loop2))
           ;;(sdiff ((MEXPT SIMP) $X 2) $X)
           ((equal (setq exp (sdiff exp (cadr z))) 0) (return 0))))
     (setq count (1- count))
     (go loop1)
     noun (return (diff%deriv (cons exp (cdr e))))))

diff

file:comm.lisp
変数の値を変えて、deriv に引数全部を渡す。
maiximaで以下を実行

diff(x^2,x,1);

;;args <- '((($F . SIMP) $X) $X 1)
(defmfun $diff (&rest args)
   ;;dynamic-extent <- '((($F . SIMP) $X) $X 1)
  (declare (dynamic-extent args))
   ;;derivlist <- nil
  (let (derivlist)
    (deriv args)))

loop

loop for 変数 in lst を使用するときはループキーワードとセットで使用する。
他の繰り返しはここを参照。
ループキーワード: append, appending, collect, collecting, nconc, nconcing はリストの中に値を蓄積し、最後にそのリストを返す。

ループキーワード: count, counting, maximize, maximizing, minimize, minimizing, sum summing は数値の累積/積算をし、最後にその値を返す。

(loop for alphabet in '(a b c d e)
      collect alphabet
      )

>(A B C D E)

;; lst に値が蓄積されるが collect はに値を返さない
(loop for alphabet in '(a b c d e)
      collect alphabet into lst
      )

>NIL

(loop for num in '(1 2 4 5 7 9 3 5 7 1)
      maximize num
      )

>9

;; append nconc の引数はlist. atomを渡すとエラー。
(loop for lst in '(a b (c d) (e f))
      if(listp lst) nconc lst)
>(C D E F)

declare-top

file:lmdcls.lisp

(defmacro declare-top (&rest decl-specs)
  `(eval-when
    ,(cond (*macro-file*  #+gcl '(compile eval load)
     #-gcl '(:compile-toplevel :load-toplevel :execute) )
    (t #+gcl '(eval compile) #-gcl '(:compile-toplevel :execute)))
    ,@(loop for v in decl-specs
      unless (member (car v) '(special unspecial)) nconc nil
      else
      when (eql (car v) 'unspecial)
      collect `(progn
         ,@(loop for w in (cdr v)
    collect #-(or gcl scl cmu ecl)
                                       `(remprop ',w
       #-excl 'special
       #+excl 'excl::.globally-special.)
    #+(or gcl scl cmu ecl)
           `(make-unspecial ',w)))
      else collect `(proclaim ',v))))

defmode

mrgmac.lisp

(declare-top (special name bas selector))

(defmacro defmode (&rest x)
  (push 'defmode x)
  (let ((selector (member 'selector (cddddr x) :test #'eq)))
    (define-mode (second x) (fourth x))
    (mapc 'eval (cddddr x))
    `',(second x)))

context

db.lisp

(defmode context ()
  (atom (selector cmark fixnum 0) (selector subc) (selector data)))

with-new-context

;; Used to temporarily bind contexts in such a way as to not cause
;; the context garbage collector to run. Used when you don't want to
;; stash away contexts for later use, but simply want to run a piece
;; of code in a new context which will be destroyed when the code finishes.
;; Note that this code COULD use an unwind-protect to be safe but since
;; it will not cause out and out errors we leave it out.

(defmacro with-new-context (sub-context &rest forms)
  `(let ((context (context ,@sub-context)))
     (prog1 ,@forms
       (context-unwinder))))