Subject: turning a sequence of epxressions into a function
From: Leo Butler
Date: Sun, 23 May 2010 17:13:26 +0100 (BST)
On Sun, 23 May 2010, Michael S. Sadler wrote:
< I am new to Maxima, and I am trying to correct the subtle things that I am
< doing wrong.
<
< If you have the time, perhaps you could examine some code for me.
<
< I have the start of a program that works, when I try to turn part of the code
< into a function, I fail.
<
< The first snippet works just fine, I am plotting a graph knowing only the
< curvature, an
< initial point and direction. The following is somewhat contrived, but it is a
< start, and works fine.
<
< I start off by defining the curvature k(a), the arclength, and a function
< called the witch of agnesi.
< I then create an array called points, big enough to hold all my increments,
< calculate the initial point and direction,
< then fill in the rest of the points using the curvature, moving along the arc
< length.
<
< this works, but I want to turn it into a function.
Welcome to Maxima.
The error you are making is a subtle one: the double quote operator ''
is implemented by the parser, so it does not work as you expect in
functions. Since you are using it to do substitution, there is an easy
fix: replace ''foo with subst('x=x,foo) in your function definition.
When I do this, I get
get(c,a,b,increment) :=
block( [len,points,f,x,y,angle],
len : arclength(c,a,b)/increment,
points:makelist([0,0],x,0,floor(len)),
f : k(c),
angle : float(atan(diff(c,x))),
x : a,
y : subst('x=x,c),
angle : subst('x=x,angle),
points[1] : [x,y],
for i : 2 thru len+1 do (
angle : float(angle),
angle : angle + float(subst('x=x,f))*increment,
x : x + cos(angle)*increment,
y : y + sin(angle)*increment,
points[i] : [x,y]
),
points
)$
and this works for me.
(Note that length is a function, and while you can create a variable
of the same name, I have changed its name).
Incidentally, dispfun(foo) will show you the Maxima definition of the
function foo. This can be useful for problems like this.
Leo
<
<
< ********************************************
<
< kill(all)$
<
< k(a) := diff(a,x,2)/(1+diff(a,x,1)^2)^(3/2)$
< arclength(f,a,b) := quad_qags(sqrt(1+diff(f,x,1)^2),x,a,b)[1]$
< witch(a) := (8*(a)^3)/(x^2+4*(a)^2)$
<
<
< c : witch(1)$
< a : -15$
< b : 15$
< increment : 0.01$
<
<
< length : arclength(c,a,b)/increment$
< points:makelist([0,0],x,0,floor(length))$
<
<
< f : k(c)$
<
< angle : float(atan(diff(c,x)))$
< x : a$
< y : ''c$
< angle : ''angle$
< points[1] : [x,y]$
<
< for i : 2 thru length+1 do (
< angle : angle + float(''f)*increment,
< x : x + cos(angle)*increment,
< y : y + sin(angle)*increment,
< points[i] : [x,y]
< )$
<
< plot2d([discrete, points]);
<
< ********************************************
<
<
< The problem that I am having is turning this into a function. Here is what I
< believe might be the right thing.
<
< ********************************************
<
< kill(all)$
<
< k(a) := diff(a,x,2)/(1+diff(a,x,1)^2)^(3/2)$
< arclength(f,a,b) := quad_qags(sqrt(1+diff(f,x,1)^2),x,a,b)[1]$
< witch(a) := (8*(a)^3)/(x^2+4*(a)^2)$
<
< get(c,a,b,increment) :=
< block( [length,points,f,x,y,angle],
< length : arclength(c,a,b)/increment,
< points:makelist([0,0],x,0,floor(length)),
< f : k(c),
< angle : float(atan(diff(c,x))),
< x : a,
< y : ''c,
< angle : ''angle,
< points[1] : [x,y],
< for i : 2 thru length+1 do (
< angle : angle + float(''f)*increment,
< x : x + cos(angle)*increment,
< y : y + sin(angle)*increment,
< points[i] : [x,y]
< ),
< points
< )$
<
< c : witch(1)$
< a : -15$
< b : 15$
< increment : 0.01$
<
< P : get(c,a,b,increment)$
<
< plot2d([discrete, P]);
<
< ********************************************
<
<
< Maxima does not like this at all. Any help would be appreciated. I think it
< may have to do with evaluating things wrong.
< _______________________________________________
< Maxima mailing list
< Maxima at math.utexas.edu
< http://www.math.utexas.edu/mailman/listinfo/maxima
<
<
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.