04/23/2009, 07:04 PM

I know you use maple, so I've tried to make this Maple-ish, but I do not have Maple, so I'm programming blind. Let me know if this works:

Andrew Robbins

Code:

`carleman_matrix := proc(expr, x, p, n) `

[seq([seq(subs(diff((expr-p)^j, x, k), x=p), k=0..n)], j=0..n)];

end proc;

iterate := proc(m, p, t)

c := m[1][1];

n := length(m);

ev := eigenvects(m);

sm := [seq(ev[3*j + 1], j=0..n)];

dt := [seq([seq(if j==k then c^k else 0 endif, k=0..n)], j=0..n)];

mt := evalm(inverse(sm) &* dt &* sm);

p + sum(mt[1][k]*(1 - p)^k, k=1..n);

end proc;

# if this works, you can use it like this:

a := 1.0001; # the base

p := 1.000999; # the fixed point

mat := carleman_matrix(a^x, x, p, 20);

flo := iterate(mat, p, 0.5);

Andrew Robbins