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 `x`

: `x0`

, `x0+h`

, and `x0-h`

. Where `h`

is 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),

'at('diff(f(x),x), x=x0),

'at('diff(f(x),x,2), x=x0),

'at('diff(f(x),x,3), x=x0),

'at('diff(f(x),x,4), x=x0),

'at('diff(f(x),x,5), x=x0)

];

f_coefs : submatrix(factor(augcoefmatrix([f_taylor], rhs_list)), 7);

fp_coefs : submatrix(factor(augcoefmatrix([fp_taylor], rhs_list)), 7);

Using the

`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)));

The matrix

`A`

represents 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

`A`

rather 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([0],[1],[0],[0],[0],[1]);

/* 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),

at(f(x),x=x0),

at(f(x),x=x0+h),

'at('diff(f(x),x),x=x0-h),

'at('diff(f(x),x),x=x0),

'at('diff(f(x),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=x0`

for 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.

The

`h**4`

term 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.

## No comments:

## Post a Comment