kill(all)$

maxima:appendix heat:special function ,version 1.0
peter.vlasschaert@gmail.com,20/07/2021

load packages,option
=================

load(contrib_ode)$
load(orthopoly)$
ratprint:false;

A\ error function (complimentary version)
================================

→ first : step 'definition error function '

myerf(x):=2/(%pi)^(1/2)*integrate(%e^(-t^2),t,0,x);

test : 

myerf(0);

myerf(∞);

assume (a>0);

myerf(-a);

→ build in (Maxima);

erf(0);

erf(∞);

erf(-a);

→ second :step 'definition (c)omplimentary  error function'

→ the same both (pa=pb)

pa:myerfc(x):=1-myerf(x);

pb:my1erfc(x):=2/(%pi)^(1/2)*integrate(%e^(-t^2),t,x,∞);

myerfc(0);

my1erfc(0);

→ build in (Maxima);

erfc(0);

→ third : step differentation.

diff(erf(x),x);

diff(erf(x),x,2);

f(x):=erf(x);

check: 'see above'

is(f(0)=erf(0));

for k:0.1 thru 5  step 1 do(
print("value = ",k)=disp(erfc(k))
    );

representation : erf 'with taylor series',big O (x^10)

taylor(erf(x),x,0,10);

B\ Bessel  Functions 
================
use : cylindrical and spherical geometries, 'heat conduction'

ref :https://maxima.sourceforge.io/docs/manual/maxima_79.html

declare(n,integer);

p1:x^2*'diff(y,x,2)+x*'diff(y,x)+(x^2-n^2)*y =0;

→ Bessel's equations of order n = integer

p1a:ode2(p1,y,x);

p2:x^2*'diff(y,x,2)+x*'diff(y,x)-(x^2+n^2)*y=0;

→modified Bessel's equations of order n

p2a:ode2(p2/x^2,y,x);

use :load packages

p2b:contrib_ode(p2,y,x);

rem :contrib_ode
use : bessel_y and bessel_i 

??  bessel_k$

 -- Function: bessel_i (<v><, <z><)
     The modified Bessel function of the first kind of order v and
     argument z.
     'bessel_i' is defined as
                              inf
                              ====   - v - 2 k  v + 2 k
                              \     2          z
                               >    -------------------
                              /     k! gamma(v + k + 1)
                              ====
                              k = 0
     although the infinite series is not used for computations.

p3b:y = %k1*bessel_i(n,x)+%k2*bessel_k(n,x);

math : 
y=C1*I[n](x)+C2*K[n](x)  solution :p2
------------------------------------------------------------------
I[n](x) , modified Bessel equation first kind.
K[n](x):modified Bessel equation second kind.
-----------------------------------------------------------------

plot function :
==========

wxplot2d([bessel_j(0,x)], [x,-5,5])$

makelist(bessel_j(k,x),k,0,5);
wxplot2d(%, [x,-5,5])$

rem : bessel_x: to plot  x:{i,jk,y}

orthopoly:
=======

15.2 Bessel Functions

Function: bessel_j (v, z)

    The Bessel function of the first kind of order v and argument z.

    bessel_j is defined as

                    inf
                    ====       k  - v - 2 k  v + 2 k
                    \     (- 1)  2          z
                     >    --------------------------
                    /        k! gamma(v + k + 1)
                    ====
                    k = 0

    although the infinite series is not used for computations.

    Categories:  Bessel functions Special functions

Function: bessel_y (v, z)

    The Bessel function of the second kind of order v and argument z.

    bessel_y is defined as

                  cos(%pi v) bessel_j(v, z) - bessel_j(-v, z)
                  -------------------------------------------
                                 sin(%pi v)

    when v is not an integer. When v is an integer n, the limit as v approaches n is taken.

    Categories:  Bessel functions Special functions

Function: bessel_i (v, z)

    The modified Bessel function of the first kind of order v and argument z.

    bessel_i is defined as

                        inf
                        ====   - v - 2 k  v + 2 k
                        \     2          z
                         >    -------------------
                        /     k! gamma(v + k + 1)
                        ====
                        k = 0

    although the infinite series is not used for computations.

    Categories:  Bessel functions Special functions

Function:  Category: Exponential Integrals

expintegral_chi expintegral_ci expintegral_e expintegral_e1 expintegral_e_simplify expintegral_ei expintegral_li expintegral_shi expintegral_si expintexpand expintrep bessel_k (v, z)

    The modified Bessel function of the second kind of order v and argument z.

    bessel_k is defined as

               %pi csc(%pi v) (bessel_i(-v, z) - bessel_i(v, z))
               -------------------------------------------------
                                      2

C\ The exponential integral function
===========================

 Category: Exponential Integrals

expintegral_chi expintegral_ci expintegral_e expintegral_e1 expintegral_e_simplify expintegral_ei expintegral_li expintegral_shi expintegral_si
expintexpand expintrep 
    when v is not an integer. If v is an integer n, then the limit as v approaches n is taken.

myEi(x):=integrate(exp(-t)/t,t,x,∞);

check values

myEi(∞);

myEi(0.001);

myEi(∞) this integral is divergent = ∞
rem:part  D , gamma function.

→ example

assume(b>0);

p3:diff(myEi(b),b);

for k:0.1 thru 5 step 1 do(
    print(k,myEi(k))
)$

D\Gamma and the Incomplete Gamma Function.
=====================================

assume(xx>0,aa>0);

→ definition (=Gamma Function), Γ(xx)

p4:integrate(exp(-t)*t^(xx-1),t,0,∞) ='integrate(exp(-t)*t^(xx-1),t,0,∞);

→ definition (=Incomplete Gamma Function), Γ(xx)

p4a:integrate(exp(-t)*t^(aa-1),t,xx,∞)='integrate(exp(-t)*t^(aa-1),t,xx,∞) ;

example :

p4b1:gamma(4);

p4b2:gamma(n+1);

p4b3::gamma(n+1);

convert : 'gamma,factorial'

makefact(gamma(n+1));

makegamma(n!);

p4b4:gamma_incomplete(1,-2);

p4b4:ev(p4b4,nouns);

p4c:gamma(1/2);

p4d:gamma(x)*gamma(1-x);

ratsimp(gamma(1-xx)*gamma(xx));

p4d = π/sin(π*x)

E\The BETA Function.
=====================================

ref : https://maxima.sourceforge.io/docs/manual/maxima_81.html

→ first :definition

mybeta(z1,z2):= (gamma(z1)*gamma(z2))/gamma(z1+z2);

mybeta(0.3,0.5);

→ second :definition

my1beta(z1,z2):= integrate((1-t)^(z1-1)*t^(z2-1),t,0,1);

float(my1beta(0.3,0.5));

build in function :beta

beta(0.3,0.5);

example :

beta_expand:true; /* properties working on the beta functions */
[beta(z1+3,z2-2),beta(z1-3,z2+2)];

rem : gamma,gamma_incomplete used in general

F\The Generalized Hypergeometric Function.
=====================================
pochhammer symbol: 'see google'

Function: hypergeometric ([a1, ..., ap],[b1, ... ,bq], x)

    The hypergeometric function. Unlike Maxima’s %f hypergeometric function, the function hypergeometric is a simplifying function; also, hypergeometric supports complex double and big floating point evaluation. For the Gauss hypergeometric function, that is p = 2 and q = 1, floating point evaluation outside the unit circle is supported, but in general, it is not supported.

    When the option variable expand_hypergeometric is true (default is false) and one of the arguments a1 through ap is a negative integer (a polynomial case), hypergeometric returns an expanded polynomial.

    Examples: 

p6:x*(1-x)*'diff(y,x,2)+(cx-(ax+bx+1)*x)*'diff(y,x)-ax*bx*y =0;

example:

p6a:hypergeometric([],[],x);

expand_hypergeometric:true;

p6b:hypergeometric([ax,bx],[cx],x);

hypergeometric([-3],[7],x);

hypergeometric([-3,4],[7],x);

hypergeometric([v,-2],[-7],x);

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