kill(all)$

maxima : taylor two variable (symbolic or discrete),version 1.0
peter.vlasschaert@gmail.com,25/02/2022

ref:https://www.youtube.com/watch?v=9H0wlv0lQl8&t=394s

load packages
-----------------

load(lrats)$

1e) integer partitions of the partial differentials 
===================================

setpartdif(m1):=block([],
/* m1: integer : 0,1,2,.......**/
load(nset),    
if m1 = 0 then s3:[0,0] 
    else
(s1:integer_partitions(m1,2),
s2:listify(s1),  
s3:copy(s2), 
     for i:1 thru length(s1) do(
          if s2[i][1]≠ s2[i][2] then
         s3:endcons(reverse(s2[i]),s3)
                                             )
  ),   
    return(s3)
)$

p1:setpartdif(2);

2e) Taylor expansion around (0,0).
==========================

a) concept:at
--------------

→ one variable

at(diff(g(x),x,2),x=a);

→ two variable

at(at('diff(f(x,y),x,1,y,1),x=x0),y=x0);

b) taylor expansion f(x,y)
-------------------------------

use : taylor from maxima

f(x,y)=taylor(f(x,y),[x,y],[x0,x0],2);

use : my taylor for maxima (around:(0,0))

yy(n,m):=ev(1/(n!*m!)*at(at(diff(f(x,y),x,n,y,m),x=x0),y=y0),x0=0,y0=0);

f(x,y)=ev(sum(sum(yy(n,m)*(x-x0)^n*(y-y0)^m,n,0,inf),m,0,inf),x0=0,y0=0);

1) first term:
=========

q1:yy(0,0);

q1aa:q1;

2) second term:
===========

q2:setpartdif(1);

q2a:makelist(yy(q2[i][1],q2[i][2])*x^q2[i][1]*y^q2[i][2],i,1,length(q2));

q2aa:sum(q2a[i],i,1,length(q2a));

3) thrid term:
===========

q3:setpartdif(2);

q3a:makelist(yy(q3[i][1],q3[i][2])*x^q3[i][1]*y^q3[i][2],i,1,length(q3));

q3a1:makelist(yy(setpartdif(2)[i][1],setpartdif(2)[i][2])*
                       x^setpartdif(2)[i][1]*y^setpartdif(2)[i][2],i,
                       1,length(setpartdif(2)));

q3aa:sum(q3a[i],i,1,length(q3a));

taylor :approximation :f(x,y) @ (0,0), until second-order derivatives.
===================================================

f(x,y)=q1aa+q2aa+q3aa;

use : my taylor for maxima (around:(x0,y0))

xx(n,m):=ev(1/(n!*m!)*at(at(diff(f(x,y),x,n,y,m),x=x0),y=y0));

f(x,y)=ev(sum(sum(xx(n,m)*(x-x0)^n*(y-y0)^m,n,0,inf),m,0,inf));

taylor :approximation :f(x,y) @ (x0,y0), until second-order derivatives.
===================================================

qq:[f(x0,y0)];

for i1:1 thru 2 do(
v[i1]:makelist(xx(setpartdif(i1)[i][1],setpartdif(i1)[i][2])*
                       (x-x0)^setpartdif(i1)[i][1]*(y-y0)^setpartdif(i1)[i][2],i,
                       1,length(setpartdif(i1))),
qq:endcons(v[i1],qq)
);

qq;

qq:flatten(qq);/* remove intern [] */

qqq:f(x,y)=sum(qq[i2],i2,1,length(qq));

set :x-x0=Δx,y-y0=Δy

qa:f(x0+Δx,y0+Δy)=expand(lratsubst([x-x0=Δx,y-y0=Δy],rhs(qqq)));

set :x=x0,y=y0

qb:subst([x0=x,y0=y],qa);

qc:hipow(qb,Δx);

qc1:length(rhs(qb));

qd:makelist(part(rhs(qb),i3),i3,1,length(rhs(qb)));

qd1a:args(qd[1]);

qe1:makelist(setpartdif(ix),ix,0,2);

qe2:flatten(qe1);

qe2:makelist([qe2[2*id-1],qe2[2*id]],id,1,6);

qe3:makelist(if qe2[iy]=0 then qe2[iy]:ss else qe2[iy]:qe2[iy],iy,1,length(qe2));

qe4:makelist([Δx,Δy],i,1,qc1);

we need to find each coef.
--------------------------------

qe5:makelist(qe4[i][1]^qe3[i][1]*qe4[i][2]^qe3[i][2],i,1,qc1);

discrete version :taylor expansion two variables
===================================

step 1: principal:
---------------------

declare_index_properties(A1,[postsubscript,postsuperscript,presubscript,presuperscript])$
A1[postsubscript,postsuperscript,presubscript,presuperscript];

declare_index_properties(A2,[postsubscript,postsuperscript])$
A2[postsubscript,postsuperscript];

A2[[x,y],[0,0]];

[x,y] = derivative x,y
[0,0] = function(A2)

A2[[x,y],[0,2]];

[x,y] = derivative x,y
[a,b] = a derivative of x, b derivative of y
a: pos integer (example: a=0)
b: pos integer (example: b=2)

step 2: makelist discrete symbols:
---------------------------------------

see above : ds1,ds2

ds1:qe3;

ds1a:part(ds1,5);

ds1b1:part(ds1,5)[1]!;
ds1b2:part(ds1,5)[2]!;

coef: for differential

ds1b:1/(ds1b1*ds1b2)!;

ds2:qe5;

declare_index_properties(yi,[postsubscript,postsuperscript])$

lhs : result

l1:yi[[0,0],[i,j]];

rhs:result

dr1:makelist(yi[[x,y],part(ds1,i)]*ds2[i],i,1,qc1);

vv:[];

for i:1 thru 6 do(
vv:endcons(yi[[x,y],part(ds1,i)]*ds2[i],vv)
);

vv;

→ correction coef ! vv or dr1

with use : ds1

ds3:makelist(1/(part(ds1,i)[1]!*part(ds1,i)[2]!),i,1,qc1);

step 3: final result (lhs=rhs)
---------------------------------------

f1:makelist(dr1[i]*ds3[i],i,1,qc1);

→ lhs=rhs

f2:l1=sum(f1[i],i,1,qc1);