First we'll use the
taylor()function to define truncated Taylor series expressions for the function and it's first derivative. Then we'll us the
ev()function to evaluate the expressions with different values of
his the step size.
f_taylor : taylor(f(x),x,x0,5);
fp_taylor : taylor('diff(f(x),x),x,x0,4);
rhs_list : [at(f(x),x=x0),
f_coefs : submatrix(factor(augcoefmatrix([f_taylor], rhs_list)), 7);
fp_coefs : submatrix(factor(augcoefmatrix([fp_taylor], rhs_list)), 7);
augcoefmatrix()function is an easy way to extract the coefficients of the various derivatives in the Taylor series. Coefficients of the function expansion:
Coefficients in the first derivative expansion:
Now we need to pack the coefficients for each point into a matrix.
/* add the columns for the points with function values: */
A : addcol(transpose(ev(f_coefs, x=x0-h)));
A : addcol(A, transpose(ev(f_coefs, x=x0)));
A : addcol(A, transpose(ev(f_coefs, x=x0+h)));
/* add the columns for the points with function derivative values: */
A : addcol(A, transpose(ev(fp_coefs, x=x0-h)));
A : addcol(A, transpose(ev(fp_coefs, x=x0)));
A : addcol(A, transpose(ev(fp_coefs, x=x0+h)));
Arepresents the linear system we need to solve to find the coefficients of the function value and its derivative to use in our finite difference approximation.
Since we might want to solve for more than one expression, it's a good idea to factor
Arather than just inverting it. The matrix should look familiar if you've used a Taylor Table before.
/* we need to invert A to get finite difference formulas so factor it: */
M : lu_factor(A);
[P,L,U] : get_lu_factors(M);
Now it is a straight forward matter to generate a right-hand-side vector for our linear system and solve for the unknown coefficients. The numerical rhs-vector represents the derivatives of the function.
Putting a one in the appropriate position of the vector gives the corresponding finite difference expression for that derivative.
/* right hand side vector for getting the first derivative: */
b_fp : matrix(,,,,,);
/* coefficients for the compact difference formula: */
a_fp : fullratsimp(lu_backsub(M, b_fp));
/* values vector (function and derivative at the points): */
vv : matrix([at(f(x),x=x0-h),
/* compact difference formula for the first derivative: */
cd_fp : vv.a_fp;
cd_fp_reduced : expand((h**4) * (cd_fp) / 30);
You might notice the extra 1 in the 6th place of the right-hand-side vector, it's there so we don't get the trivial solution of using the value for the derivative at
x=x0for our estimate of the value of the derivative (we obviously don't know it yet). This results in the following difference expression.
Which can be multiplied by an arbitrary constant to achieve a more computation ready form.
h**4term can be neglected, resulting in the well-known fourth order compact (implicit) expression for the first derivative.
Here's the Maxima batch file used for this.