axiom-mail
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Axiom-mail] A tiny type conversion problem


From: Martin Rubey
Subject: Re: [Axiom-mail] A tiny type conversion problem
Date: 19 Jul 2008 09:18:26 +0200
User-agent: Gnus/5.09 (Gnus v5.9.0) Emacs/21.4

"Alasdair McAndrew" <address@hidden> writes:

> I'm a bit embarrassed not to be able to work this out for myself, but,
> anyway...  Now that I've discovered the joys of the NumericalQuadrature
> package, I want to use it to calculate the total variation of some
> functions; the total variation being defined (for differentiable functions)
> as the integral of the absolute value of the derivative.  So I'd like to use
> something like
> 
> romberg(x+->abs(D(x/(x^3+1),x)),0.0,1.0,0.1,0.1,10,10)
> 
> But this doesn't work.  Romberg requires an anonymous function float ->
> float, but the use of D returns an expression. (Expression integer, in this
> example).  How, and where, do I place the type conversions so that romberg
> will do its job?

Although it looks ugly, it is in fact quite mathematical: you have to evaluate
the function at the givel point.  You can either use quotes to indicate the
places where you mean the symbol x rather than the place:

(7) -> romberg(x+->abs eval(D('x/('x^3+1),'x), 'x=x),0.0,1.0,0.1,0.1,10,10)

   (7)
   [value= 0.5582671916 6951839528,error= 0.0,totalpts= 1025,success= false]
   Type: Record(value: Float,error: Float,totalpts: Integer,success: Boolean)

or simply use a different name:

(8) -> romberg(x+->abs eval(D(y/(y^3+1),y), y=x),0.0,1.0,0.1,0.1,10,10)

   (8)
   [value= 0.5582671916 6951839528,error= 0.0,totalpts= 1025,success= false]
   Type: Record(value: Float,error: Float,totalpts: Integer,success: Boolean)

actually, you want success:

(9) -> romberg(x+->abs eval(D(y/(y^3+1),y), y=x),0.0,1.0,0.1,0.1,10,20)

   (9)
   [value= 0.5582674141 9162960882, error= 0.2225221112 135 E -6,
    totalpts= 2049, success= true]
   Type: Record(value: Float,error: Float,totalpts: Integer,success: Boolean)

Let's check:

(10) -> (integrate(D(y/(y^3+1),y),y=0..(1/2)^(1/3), 
"noPole")-integrate(D(y/(y^3+1),y),y=(1/2)^(1/3)..1, "noPole"))::EXPR Float

   (10)  0.5582673679 7879964984
                                                       Type: Expression Float

Martin





reply via email to

[Prev in Thread] Current Thread [Next in Thread]