kill(all)$

calculus : Dirac Delta of function , version 1.0
peter.vlasschaert@gmail.com,11/02/2020

load "packages"
============

load(lrats)$

integral representation : Delta(x-c)
------------------------------------------

a)  range : -∞,∞
--------------------

p1a:integrate(delta(x-c)*f(x),x,minf,inf) =f(c);

b) range :-c-ε,+c+ε
-------------------------

c >>> ε

p1b:integrate(delta(x-c)*f(x),x,-c-ε,+c+ε) =f(c);

find : delta(g(x))
============
use : g , continues differentiable function with g' nowhere zero

p2:g(x) = taylor(g(x), x, a, 1);

p2a:lhs(p2) =rhs(p2),simpsum;

p2b:lhs(p2) = factorsum(rhs(p2a));

because:   a: 'root of g(x)' → g(a) =0
use :truncated taylor series , @a
=>  a  must be real root from, g(x)

p2c:subst([g(a)= 0],p2b);

p2c1:rhs(p2c) = lhs(p2c);

rem : constraint :g'(a) ≠ 0

use : princepal p1a

p3:'integrate(f(x)*delta(g(x)),x,minf,inf);

p3a:subst([p2c],p3);

p4:p3=p3a  ;

p4a:args(lhs(p4))  = args(rhs(p4));

p4b:lhs(p4a)[1] =rhs(p4a)[1];

p4c1:c=at('diff(g(x),x,1),x=a);

p4c2:rhs(p4c1) = lhs(p4c1);

c = value , constant 

p4c3:subst([p4c2],p4b)/f(x);

we want rule :
cc = constant

scaling property : dirac delta function

rule : scaling property 'dirac delta'

p5:delta(cc*x) =1/abs(cc)*delta(x);

p4c4:subst([cc=c, x=x-a],p5);

final result : delta(g(x)), g(a) = 0, only one root: x=a
======================================

p4c5:subst([p4c1],p4c4);

p4c6:subst([p2c1],lhs(p4c5)) = rhs(p4c5);

example : delta(k*x+b)
================

p6:delta(k*x+b) =delta(k*m);

p6a:args(lhs(p6))[1] =args(rhs(p6))[1];

p6b:expand(solve(p6a,m))[1];

constraint : k ≠ 0

use :p5

p6c:subst([cc=k,x=rhs(p6b)],p5);

example : delta(x^2-a^2)
==================

p7:g(x) = x^2-a^2;
p7a:lhs(p7)=factor(rhs(p7));

1e) find roots
==========

p7b:solve(rhs(p7),x);

2e) differentation
=============

p7c:diff(p7a,x);

root : x=a
=======

p7c1:at(p7c,x=a);

p7c1a:abs(p7c1);

root : x=-a
========

p7c2:at(p7c,x=-a);

p7c2a:abs(p7c2);

3e) integral :sum, superposition
=======================

roots :x=a

p7c3a:part(rhs(p7a),1);

p7c3a1:delta(p7c3a);

roots : x=-a

p7c3b:part(rhs(p7a),2);

p7c3b1:delta(p7c3b);

4e) final result
===========

p7c4a:delta(x^2-a^2)= 1/rhs(p7c1a)*p7c3a1+1/rhs(p7c2a)*p7c3b1;

p7c4b:lhs(p7c4a)=combine(rhs(p7c4a));

when : g(x), a[i] real roots
===================
see example : above. delta(x^2-a^2)., a[1]=a,a[2]=-a.

p8:lhs(p4c6)=sum(subst([a=a[j]],rhs(p4c6)),j,1,i);

p8a:args(p8);

p8a1:p8a[2];

p8a2:args(p8a1);

p8a21:denom(p8a2[1]);

p9a:integrate(f(x)*delta(g(x)),x,minf,inf);

p9b:sum(f(a[j])/p8a21,j,1,i);

integral :result

p9:p9a=p9b;

example : i=2, x^2-a^2;a[1] = a,a[2] = -a
-------------------------------------------------

 ;

p10:subst([i=2,g(x)=x^2-a^2],p9);

p10a:p10,simpsum;

p10b:ev(p10a,nouns);

p10b1:subst([a[1]=a,a[2]=-a],p10b);

p10b2:lhs(p10) = combine(rhs(p10b1));