On 3/30/07, Zoho Vignochi <zoho.vignochi at gmail.com> wrote:
> believe the problem is that "step" returns 2 values and I am assigning
> that to ans and calling each value as ans[1] and ans[2]. Perhaps that is
> not legal in maxima.
It is legal to return [a,b] and reference that value as val[1] and val[2].
I have not looked at your program in detail, but there are five simple
errors you'll need to correct....
You use ":=" as the assignment statement e.g. a0[i]:=... This should
be just ":"; in Maxima, ":=" is used for defining functions.
Secondly, in a couple of places, you write
for I:1 thru N-1 do
this,
that,
theother,
or
while ... do
this,
that,
theother,
apparently expecting the indentation to determine the scope of the
"do". But Maxima ignores indentation; the scope of the "do" is
exactly one statement. I think what you mean is
for I:1 thru N-1 do
( this,
that,
theother ),
...
You can't use "step" as a function name, as it is a keyword in
constructions like (for i STEP j do...).
Finally, you appear to be using
diag_points:[(diag_points),[lam,alpha]]
to accumulate a list of diag_points. But that will create a nested
list. I think you want
diag_points: cons([lam,alpha], diag_points)
Also, you appear to have missed converting ap[3][0] to ap3[0] in a
couple of places.
I'm including my modified code below. I have *not* looked at its
mathematical functionality, just its Maxima syntax...
Hope this helps.
-s
-----------------------------------
n:30$
stepp(lam,alpha):=
block ([h:1/(2*n)],
a0[0]:alpha,
a1[0]:0,
a2[0]:-lam/2*sin(a0[0]),
a3[0]:0,
ap0[0]:0,
ap1[0]:0,
ap2[0]:-0.5*sin(alpha),
ap3[0]:0,
for i thru n do (
a0[i]: a0[i-1]+a1[i-1]*h+a2[i-1]*h^2,
a1[i]: a1[i-1]+2*a2[i-1]*h+3*a3[i-1]*h^2,
a2[i]: -lam/2*sin(a0[i]),
a3[i]: -lam/6*cos(a0[i])*a1[i],
ap0[i]: ap0[i-1]+ap1[i-1]*h+ap2[i-1]*h^2,
ap1[i]: ap1[i-1]+2*ap2[i-1]*h+3*ap3[i-1]*h^2,
ap2[i]: -0.5*sin(a0[i])-lam/2*cos(a0[i])*ap0[i],
ap3[i]:
-1/6*cos(a0[i])*a1[i]+lam/6*sin(a0[i])*ap0[i]*a1[i]-lam/6*cos(a0[i])*ap1[i]),
[ a0[n], ap0[n] ] ) $
newton(lam_guess,alpha):=
block ( [lam, lam_old, ans],
lam:lam_guess,
lam_old:lam+1,
while abs(lam-lam_old) > 10^(-8) do (
ans: stepp(lam,alpha),
lam_old: lam,
lam: lam-ans[1]/ans[2]),
lam)$
diagram(m) :=
block ( [lam, alpha, N:30],
lam:float(m^2*%pi^2),
diag_points:[[lam,0]],
for i thru N-1 do (
alpha:float(i/N*%pi),
lam:newton(lam,alpha),
diag_points: cons([lam,alpha], diag_points) ),
plot2d(diag_points))$