kill(all)$

secondorder finite difference (approximation) , version 1.0
peter.vlasschaert@gmail.com,06/04/2019

How to approximate  second order differentation.
=========================================

p1: 'diff(f,x,2) ;

taylor expansion around  'x+x0'  , general n -> infinity
============================================

p2:f(x+x0) =sum('diff(f(x),x,i)/i!*(x0)^i,i,0,n);

we only need  '4' terms for approximation, taylor expansion.
==================================================

p2a:subst([n = 3], p2);

we can't use here : function discriptor or define.
------------------------------------------------------------------------
reason : f() not possible on both sides.

1e) how to expand : p2a, four terms
===========================
use : command : 'simpsum'

p2b:p2a,simpsum;

2e) x0 replaced by Δx
==================

p2c:ev(p2b,x0=Δx);

approximation 'use p2c'  for  'p1'
===========================

1e)  f(x+Δx)   ?
============

p3a: p2c;

2e)  f(x-Δx)   ?
============

p3b:subst([Δx=-Δx],p2c);

3e)  1e)  + 2e)
===========

p3c:p3a+p3b;

4e) solve for 'diff(f(x),x,2)
=====================

p3d:solve(p3c,'diff(f(x),x,2))[1];

we don't have error term   'how to find error term'
=========================================

solution : 5 terms ,before only 4 terms
================================

1e)  n =4  into  'p2'
===============

p4a:subst([n=4],p2),simpsum;

2e) x0 by Δx
===========

p4b:subst([x0=Δx],p4a);

3e) x0 by -Δx
===========

p4c:subst([x0=-Δx],p4a);

4e)  2e)+3e)
==========

p4d:p4b+p4c;

5e) solve for 'diff(f(x),x,2)
=====================

p4e:solve(p4d,'diff(f(x),x,2))[1];

p4f:expand(p4e);

6e) error term
============

p4g:error_term=rhs(p4f)-rhs(p3d);

p4h:ratsimp(p4g);

general: approximation , second order differentation.
=============================================

p5:p1 = rhs(p3d)+lhs(p4h);

how to change representation 
=========================

p5a:f(x+Δx) = f[i+1];

p5b:f(x) = f[i];

p5c:f(x-Δx) = f[i-1];

p5d:[p5a,p5b,p5c];

function : one variable f(x)
======================

p6:subst(p5d,p5);

function : two variables  u(x,y)
=========================

p6a:subst([f=uu],p6);

 diff(u,x,2) ?
==========

change representation for two variables
-------------------------------------------------------------

p6b:[uu[i+1] =u[i+1,j],uu[i] = u[i,j],uu[i-1]=u[i-1,j]];

p6c:subst([uu=u],lhs(p6a)) = subst(p6b,rhs(p6a));

p6d:lhs(p4h) = subst([f(x)=u(x,y)],rhs(p4h));

 diff(u,y,2)?
=========

change representation for two variables
-------------------------------------------------------------

p7b:[uu[i+1] =u[i,j+1],uu[i] = u[i,j],uu[i-1]=u[i,j-1],x=y];

p7c:subst([uu=u,x=y],lhs(p6a)) =subst(p7b,rhs(p6a));

p7d:lhs(p4h) = subst([f(x)=u(xx,yy)],rhs(p4h))$

p7e:lhs(p4h)=subst([x=y],rhs(p7d))$

p7:subst([Δx=Δy,xx=x,yy=y],p7e);

info  :
====

between  (j-1,j,j+1)  = 2*Δy
between  (i-1,i,i+1)  = 2*Δx
=======================
Δx  = grid unit ,x -direction.
Δy  = grid unit , y -direction.
i index : x-direction
j index : y-direction
=======================