Hello,
Thank you very much for your help. I am still missing something however.
I took your advice and altered the program as you suggested and I am now
here
On Fri, 2007-30-03 at 11:48 -0700, Richard Fateman wrote:
> It could look like this, by ELIMINATING ALL DECLARATIONS..
> step(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), ap[3][0]:0,
>
...
> m/6*cos(a0[i])*ap1[i],
> these[i]:=[op(these),[i*h,a0[i]]], /* these is not used here*/
> [a0[n],ap0[n]])
n:30,
step(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),
ap[3][0]:0,
for i:1 step 1 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*ap[3][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:step(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:1 step 1 thru N-1 do
alpha:float(i/N*%pi),
lam:newton(lam,alpha),
diag_points:[(diag_points),[lam,alpha]],
plot2d(diag_points))
I 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. Any suggestions?