Thanks; using your fix of replacing ev() with subst() and adding f to
either the argument list (as in your example), or to the local
variable list of the block fixes things, no warnings or errors. I
guess it's a matter of preference where you declare f, I tend not to
know my f's analytically, so no need to pass it as an argument.
explicit_fde(steps, deriv) := block([A, a_f, b, coef_list, fde,
func_vals, f_coefs, f_taylor, M, n, points, i, f],
n : length(steps),
f_taylor : taylor("f"(x), x, x0, n),
points : [x=x0+steps[1]*h],
for i : 2 thru n do
(
points : append(points, [x=x0+steps[i]*h])
),
/* right-hand side vector (based on the derivative requested): */
b : zeromatrix(n,1),
b[deriv+1,1] : 1,
/* function values vector: */
func_vals : matrix([at(f(x), points[1])]),
for i : 2 thru n do
(
func_vals : addcol(func_vals, [at(f(x), points[i])])
),
/* coefficients list */
coef_list : [at(f(x),x=x0)],
for i : 2 thru n do
(
coef_list : append(coef_list, [ 'at('diff(f(x),x,i-1), x=x0) ] )
),
f_coefs : submatrix(factor(augcoefmatrix([f_taylor], coef_list)), n+1),
/* system matrix: */
A : addcol(transpose(ev(f_coefs, points[1]))),
for i : 2 thru n do
(
A : addcol(A, transpose(subst(points[i],f_coefs)))
),
M : lu_factor(A),
/* solve for the coefficients: */
a_f : fullratsimp( lu_backsub(M, b) ),
/* finite difference expression: */
fde : fullratsimp(func_vals.a_f),
return(fde)
);
2009/5/22 Martin Sch?necker <ms_usenet at gmx.de>:
> Joshua Stults schrieb:
>>
>> Hello,
>>
>> I'm trying to learn how to use compile_file(), but am having a bit of
>> trouble.
>
> I have no experience with compliling either, but let's see: Firstly, it
> seems that the compiler cannot handle the undefined function f, however
> doing as maxima proposes:
>
>> Warning: f has a function or macro call which has not been translated
>> properly.
>> The function was totaly undefined. Maybe you want to quote it.
>
> ...doesn't seem to work. ?So we try to give the function name as a third
> argument (There may be more elegant approaches?). ?Then, i does not seem to
> be evaluated to the number it is ought to have in the for loop. ?This seems
> to be due to the evaluation partialities of ev() (of which I'm unaware). ?It
> seems that the correct position of points should be substituted into
> f_coefs, so we'll try to use subst().
>
> Now the following compiles and seems to work properly (om Maxima 5.18.1):
>
> explicit_fde(steps, deriv, f) := block([A, a_f, b, coef_list, fde,
> func_vals, f_coefs, f_taylor, M, n, points, i],
> ?n : length(steps),
> ?f_taylor : taylor(f(x), x, x0, n),
> ?points : [x=x0+steps[1]*h],
> ?for i : 2 thru n do
> ?(
> ? ?points : append(points, [x=x0+steps[i]*h])
> ? ?),
> ?/* right-hand side vector (based on the derivative requested): */
> ?b : zeromatrix(n,1),
> ?b[deriv+1,1] : 1,
> ?/* function values vector: */
> ?func_vals : matrix([at(f(x), points[1])]),
> ?for i : 2 thru n do
> ?(
> ? ?func_vals : addcol(func_vals, [at(f(x), points[i])])
> ? ?),
> ?/* coefficients list */
> ?coef_list : [at(f(x),x=x0)],
> ?for i : 2 thru n do
> ?(
> ? ?coef_list : append(coef_list, [ 'at('diff(f(x),x,i-1), x=x0) ] )
> ? ?),
> ?f_coefs : submatrix(factor(augcoefmatrix([f_taylor], coef_list)), n+1),
> ?print(f_coefs),
> ?/* system matrix: */
> ?A : addcol(transpose(ev(f_coefs, points[1]))),
> ?for i : 2 thru n do
> ?(
> ? ?A : addcol(A, transpose(subst(points[i], f_coefs)))
> ? ?),
> ?M : lu_factor(A),
> ?/* solve for the coefficients: */
> ?a_f : fullratsimp( lu_backsub(M, b) ),
> ?/* finite difference expression: */
> ?fde : fullratsimp(func_vals.a_f),
> ?return(fde)
> ?);
>
>
> (%i1) kill(all);
> (%o0) done
> (%i1) compile_file("m:/maxima/explicit_fde.mac");
> Translation begun on #pm:/maxima/explicit_fde.mac.
> See the `unlisp' file for possible optimizations.
> Compiling m:/maxima/explicit_fde.LISP.
> End of Pass 1.
> End of Pass 2.
> OPTIMIZE levels: Safety=2, Space=3, Speed=3
> Finished compiling NIL.
> (%o1)
> [m:/maxima/explicit_fde.mac,m:/maxima/explicit_fde.LISP,m:/maxima/explicit_fde.UNLISP,#pm:/maxima/explicit_fde.o]
> (%i2) explicit_fde([-2, -1,0,1,2], 1, f);
> matrix([1,-(x0-x),(x0-x)^2/2,-(x0-x)^3/6,(x0-x)^4/24])
> Is ?h ?positive, negative, or zero?p;
> (%o2) -(f(x0+2*h)-8*f(x0+h)+8*f(x0-h)-f(x0-2*h))/(12*h)
>
--
Joshua Stults
Website: http://j-stults.blogspot.com