Changeset 15460 in project


Ignore:
Timestamp:
08/14/09 09:59:55 (10 years ago)
Author:
Ivan Raikov
Message:

improvements to the differentiation code in nemo-utils

Location:
release/4/nemo/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • release/4/nemo/trunk/examples/carelli05_run.m

    r15451 r15460  
    11#!/usr/bin/env octave
    22
    3 function ys = be_newton ( f, fy, t, y_initial, h )
     3function ys = beuler ( f, fy, t, y_initial, h )
    44%
    5 %  function [ x, y ] = be_newton ( f, fy, t, y_initial, h )
     5%  function ys = beuler ( f, fy, t, y_initial, h )
    66%
    7 BE_NEWTON uses NSTEP steps of the backward Euler Newton method
    8 %  to estimate Y, the solution of an ODE, at the equally spaced points X in
    9 %  the range X_RANGE(1) to X_RANGE(2). 
     7Uses the backward Euler method
     8%  to estimate Y, the solution of an ODE, at equally spaced points in
     9%  the range T(1) to T(2). 
    1010%
    1111%  Euler's forward method is used as a predictor,
  • release/4/nemo/trunk/nemo-matlab.scm

    r14732 r15460  
    672672         (pp indent ,nl (endfunction))
    673673
     674         (pp indent ,nl (function dy = ,(s+ sysname "_dy")
     675                                  (,(sl\ ", " `(t y)))))
     676
     677         
     678
     679         (pp indent ,nl (endfunction))
     680
     681
    674682         (if (not (null? steady-state-index-map))
    675683             (begin
  • release/4/nemo/trunk/nemo-utils.scm

    r15459 r15460  
    221221(define LOG2E  1.44269504088896)
    222222
    223 (define (differentiate fenv x t)
     223(define (differentiate fenv xenv x t)
    224224  (cond ((number? t)  0.0)
    225         ((symbol? t)  (if (equal? x t) 1.0 0.0))
     225        ((symbol? t)  (cond ((equal? x t) 1.0)
     226                            ((member t xenv) t)
     227                            (else 0.0)))
    226228        (else (match t
    227                 (('neg u)  `(neg ,(differentiate fenv x u)))
    228 
    229                 (('+ u v)  `(+ ,(differentiate fenv x u) ,(differentiate fenv x v)))
    230                 (('- u v)  `(- ,(differentiate fenv x u) ,(differentiate fenv x v)))
    231 
    232                 (('* (and u (? number?)) v)    `(* ,u ,(differentiate fenv x v)))
    233                 (('* v (and u (? number?)))    `(* ,u ,(differentiate fenv x v)))
    234 
    235                 (('* u v)     `(+ (* ,(differentiate fenv x u) ,v)
    236                                   (* ,u ,(differentiate fenv x v))))
    237 
    238                 (('/ u v)     `(/ (- (* ,(differentiate fenv x u) ,v)
    239                                      (* ,u ,(differentiate fenv x v)))
     229                (('neg u)  `(neg ,(differentiate fenv xenv x u)))
     230
     231                (('+ u v)  `(+ ,(differentiate fenv xenv x u) ,(differentiate fenv xenv x v)))
     232                (('- u v)  `(- ,(differentiate fenv xenv x u) ,(differentiate fenv xenv x v)))
     233
     234                (('* (and u (? number?)) v)    `(* ,u ,(differentiate fenv xenv x v)))
     235                (('* v (and u (? number?)))    `(* ,u ,(differentiate fenv xenv x v)))
     236
     237                (('* u v)     `(+ (* ,(differentiate fenv xenv x u) ,v)
     238                                  (* ,u ,(differentiate fenv xenv x v))))
     239
     240                (('/ u v)     `(/ (- (* ,(differentiate fenv xenv x u) ,v)
     241                                     (* ,u ,(differentiate fenv xenv x v)))
    240242                                  (pow ,v 2.0)))
    241243
    242                 (('cube u)     (differentiate fenv x `(pow ,u 3.0)))
    243 
    244                 (('pow u n)    (chain fenv x u `(* ,n (pow ,u (- ,n 1.0)))))
     244                (('cube u)     (differentiate fenv xenv x `(pow ,u 3.0)))
     245
     246                (('pow u n)    (chain fenv xenv x u `(* ,n (pow ,u (- ,n 1.0)))))
    245247               
    246                 (('sqrt u)     (chain fenv x u `(/ 1.0 (* 2.0 (sqrt ,u)))))
    247 
    248                 (('exp u)      (chain fenv x u `(exp ,u)))
    249 
    250                 (('log u)      (chain fenv x u `(/ 1.0 ,u)))
    251 
    252                 (('log10 u)    (chain fenv x u `(* ,LOG10E (/ ,(differentiate fenv x u) ,u))))
     248                (('sqrt u)     (chain fenv xenv x u `(/ 1.0 (* 2.0 (sqrt ,u)))))
     249
     250                (('exp u)      (chain fenv xenv x u `(exp ,u)))
     251
     252                (('log u)      (chain fenv xenv x u `(/ 1.0 ,u)))
     253
     254                (('log10 u)    (chain fenv xenv x u `(* ,LOG10E (/ ,(differentiate fenv xenv x u) ,u))))
    253255       
    254                 (('log2 u)     (chain fenv x u `(* ,LOG2E (/ ,(differentiate fenv x u) ,u))))
    255 
    256                 (('log1p u)    (differentiate fenv x `(log (+ 1.0 ,u))))
    257 
    258                 (('ldexp u n)  (differentiate fenv x `(* ,u ,(expt 2 n))))
     256                (('log2 u)     (chain fenv xenv x u `(* ,LOG2E (/ ,(differentiate fenv xenv x u) ,u))))
     257
     258                (('log1p u)    (differentiate fenv xenv x `(log (+ 1.0 ,u))))
     259
     260                (('ldexp u n)  (differentiate fenv xenv x `(* ,u ,(expt 2 n))))
    259261       
    260                 (('sin u)      (chain fenv x u `(cos ,u)))
     262                (('sin u)      (chain fenv xenv x u `(cos ,u)))
    261263               
    262                 (('cos u)      (chain fenv x u `(neg (sin ,u))))
    263 
    264                 (('tan u)      (differentiate fenv x `(* (sin ,u) (/ 1.0 (cos ,u)))))
     264                (('cos u)      (chain fenv xenv x u `(neg (sin ,u))))
     265
     266                (('tan u)      (differentiate fenv xenv x `(* (sin ,u) (/ 1.0 (cos ,u)))))
    265267               
    266                 (('asin u)     (chain fenv x u `(/ 1.0 (sqrt (- 1.0 (pow ,u 2.0))))))
     268                (('asin u)     (chain fenv xenv x u `(/ 1.0 (sqrt (- 1.0 (pow ,u 2.0))))))
    267269                                             
    268                 (('acos u)     (chain fenv x u `(/ (neg 1.0) (sqrt (- 1.0 (pow ,u 2.0))))))
     270                (('acos u)     (chain fenv xenv x u `(/ (neg 1.0) (sqrt (- 1.0 (pow ,u 2.0))))))
    269271                                             
    270                 (('atan u)     (chain fenv x u `(/ 1.0 (+ 1.0 (pow ,u 2.0)))))
     272                (('atan u)     (chain fenv xenv x u `(/ 1.0 (+ 1.0 (pow ,u 2.0)))))
    271273                                             
    272                 (('sinh u)     (differentiate fenv x `(/ (- (exp ,u) (exp (neg ,u))) 2.0)))
    273 
    274                 (('cosh u)     (differentiate fenv x `(/ (+ (exp ,u) (exp (neg ,u))) 2.0)))
     274                (('sinh u)     (differentiate fenv xenv x `(/ (- (exp ,u) (exp (neg ,u))) 2.0)))
     275
     276                (('cosh u)     (differentiate fenv xenv x `(/ (+ (exp ,u) (exp (neg ,u))) 2.0)))
    275277               
    276                 (('tanh u)     (differentiate fenv x `(/ (sinh ,u) (cosh ,u))))
    277 
    278                 (('let bnds body)  `(let ,(map (match-lambda ((v b) `(,v ,(differentiate fenv x b))))
    279                                                bnds)
    280                                       ,(differentiate fenv x body)))
    281                                                              
     278                (('tanh u)     (differentiate fenv xenv x `(/ (sinh ,u) (cosh ,u))))
     279
     280                (('let bnds body)  (let ((xenv1 (fold (lambda (vb xenv)
     281                                                        (cons (car vb) xenv)) xenv bnds)))
     282                                     `(let ,(map (match-lambda ((v b) `(,v ,(differentiate fenv xenv x b))))
     283                                                 bnds)
     284                                        ,(differentiate fenv xenv1 x body))))
    282285
    283286                ((op u)        (cond ((lookup-def op fenv) =>
    284                                       (match-lambda ((fx) (chain fenv x u `(,fx ,u)))
     287                                      (match-lambda ((fx) (chain fenv xenv x u `(,fx ,u)))
    285288                                                    (else #f)))
    286289                                     (else #f)))
     
    288291                ((op . us)     (cond ((lookup-def op fenv) =>
    289292                                      (lambda (fs)
    290                                         `(+ . ,(map (lambda (fu u) (chain fenv x u `(,fu ,u)))
     293                                        `(+ . ,(map (lambda (fu u) (chain fenv xenv x u `(,fu ,u)))
    291294                                                    fs us))))
    292295                                     (else #f)))
     
    294297                (else          #f)))))
    295298
    296 (define (chain fenv x t u)
     299(define (chain fenv xenv x t u)
    297300  (if (symbol? t) u
    298       `(* ,(differentiate fenv x t) ,u)))
     301      `(* ,(differentiate fenv xenv x t) ,u)))
    299302
    300303
     
    325328         
    326329         (('let bnds body)
    327           `(let ,(map (match-lambda ((v b) `(v ,(recur b)))
     330          `(let ,(map (match-lambda ((v b) `(v ,(simplify b)))
    328331                                    (else #f)) bnds)
    329              ,(recur body)))
     332             ,(simplify body)))
    330333         
    331334         ((op . ts) 
Note: See TracChangeset for help on using the changeset viewer.