kill(all)$

maxima : approximation 1e order differentation f(x) = f'(x) with fft,version 1.0
peter.vlasschaert@gmail.com, 13/01/2021

ref :see website : https://petervlaschemeng.weebly.com/
------------------------------------------------------------------------------------------------------------
→ SPECIAL TOPICS  :
 
=>  appendix:

-------------------------------------------------------------------------------------------------------------
new: 04 /01/2021:maxima: relation between fft & inverse_fft ,version 1.0
download  (Maxima,code)       :maximabasicfftandinverseifftversion1
download  (Maxima,png)         :maximabasicfftandinverseifftversion1outp
new: 06 /01/2021:maxima: relation between fft & inverse_fft (1D,2D),version 1.0a
download  (Maxima,code)       :maximabasicfftandinverseifftversion1a
download  (Maxima,png)         :maximabasicfftandinverseifftversion1aoutp
-------------------------------------------------------------------------------------------------------------

load packages ,option
------------------------------

load(lrats)$ load(fft)$
ratprint:false $

info :
====
fft               : fast fourier transformation
inverse_fft : inverse fast fourier transformation

  differentation : function (periodic)
===========================

pp:5;           /* value for the power of 2 */

L:2*%pi;       /* length : interval */
N: 2^pp;        /* number of points : interval */
xv[cv]:=L/N*cv; /* values : L ,interval */

function       : → f(x)

xl:makelist(xv[i],i,0,N-1)$ /* list :values */

f(xx):=sin(xx);

differentation :→ d(f(x))/dx : exact

fd(xx):=ev(diff(f(xx),xx),nouns);

result : exact differentation of f(x)

→ symbolic

---------------------------------------------

fda:fd(value);

---------------------------------------------

→ numerical ( = list of values)

--------------------------------------------

fdl:makelist(float(subst([value=xv[i]],fd(value))),i,1,N)$

--------------------------------------------

fl:makelist(f(xv[i]),i,0,N-1),float$ /* list:values of function */

length(fl); /* fl : number values */

→ fft  ( function).

f2:fft(fl)$  /* fft :fast fourier transform */

→ wave number ( = k)

p1:makelist(i,i,0,N/2-1)$
p2:[0]$
p3:makelist(i-N,i,N/2+1,N-1)$

k1:(-1)*flatten(endcons(p3,endcons(p2,p1)))$

pp=2, N=8 :  k1 = [0,-1,-2,-3,0,3,2,1]

k:((2*%pi*%i)/L)*k1 $
k:float(k)$

→ differentation : function ( = approximation )

1e)  product         : k*f2

dp:expand(makelist(f2[i]*k[i],i,1,N)),float$

2e)  differentation : function 

df:inverse_fft(dp)$ /* inverse from fft */

rdf:realpart(df)$   /* realpart :inverse_fft */

lmi:lmin(rdf); /* minimum y value : vertical*/
lma:lmax(rdf); /* maximum y value : vertical*/

info : horizontal : x values on the interval

xydaprox:realpart(makelist([xv[i],rdf[i]],i,1,length(fl)))$

3e) a)  plot               : 'differentation : function' ,approximation 

xf:['discrete,xydaprox]$

xf1:[style,lines,points], [color,blue,red]$

wxplot2d([xf,xf],xf1,[y,lmi,lma],grid2d);

b)  plot               : 'differentation : function' ,exact 

xyda:makelist([xv[i],fdl[i]],i,1,length(fl))$

wxplot2d(['discrete,xyda],[y,lmi,lma],grid2d);

combination :plot ' original function , derivative of function '

see: fl,xv

xyf:makelist([xv[i],fl[i]],i,1,length(fl))$

xp1:['discrete,xyda]$
xp2:['discrete,xyf]$
xq1:[title, "combination:f(x) & f'(x),with fft "]$
xq2:[legend, "f'(x)", "f(x)"]$
xq3:[xlabel, "x values"]$
xq4:[ylabel, "f'(x),f(x)"]$

lmif:lmin(fdl); /* minimum y value (function) : vertical*/
lmaf:lmax(fdl); /* maximum y value (function) : vertical*/

rem : combination : find : max ,min  y values

lmicomb:lmin([lmif,lmi]);
lmacomb:lmax([lmaf,lma]);

wxplot2d([xp1,xp2],[y,lmicomb,lmacomb],
    xq1,xq2,xq3,xq4,grid2d);

rem : for good plot : change the L value. 

wave numbers  : ( = k)
----------------------
→  L = interval 

maxima : k ( see : above k = 2*π/L*k1),numerical
======== 
N:8
k1 => [0,-1,-2,-3,0,3,2,1]
----------------------------------------------------------------------------------
p1:makelist(i,i,0,N/2-1)$
p2:[0]$
p3:makelist(i-N,i,N/2+1,N-1)$
k1:(-1)*flatten(endcons(p3,endcons(p2,p1)))$
----------------------------------------------------------------------------------

matlab : k ( see : above k = 2*π/L*k1),numerical
======== 
N = 8
k1 => [0,  1,  2,  3,  0,  -3,  -2,  -1]
-----------------------------------------------------------------------------------
check this with : Euler Math (=free) or Matlab
>matlab on
>N=8;
>k1 = [0:N/2-1, 0, -N/2+1:-1]
 [0,  1,  2,  3,  0,  -3,  -2,  -1]
-----------------------------------------------------------------------------------

python : k ( see : above k = 2*π/L*k1),numerical
======== 

-----------------------------------------------------------------------------------
if N % 2 == 0:
k = np.asarray(range(0, N / 2) + [0] + range(-N / 2 + 1,0))
else:
k = np.asarray(range(0,(N - 1) / 2) + [0] + range(-(N - 1) / 2, 0))
k *= 2 * np.pi / L
-----------------------------------------------------------------------------------


info:
------
maxima:

========================================================
https://maxima.sourceforge.io/docs/manual/maxima_107.html

packages :
https://maxima.sourceforge.io/docs/manual/
maxima_364.html#Category_003a-Package-fft
========================================================


matlab:

========================================================
Euler math : python mode , matlab mode ,maxima mode
ref :http://euler-math-toolbox.de/download.html
========================================================

python:

========================================================
https://kitchingroup.cheme.cmu.edu/blog/2013/02/26/Derivatives-by-FFT/
========================================================