Vector calculation of HermiteH Polynomials

The function im trying to code is hermiteH and it is as follows:

note: this is Wolfram syntax.

I want to create a vector of n hermite polynomials, and then each one to go through, some sort of matrix operation comes to mind but the implementation eludes me.

note: taken from wikipedia: Hermite polynomials.

My goal is to get something like be run through.

In the example ill get with .

My current Code:

Hk = @(k,x) hermiteH(k,x/sqrt(2))*2.^(-k/2);
k = 2;
hPoly = zeros(1,k+1);
for i=1:k
    hPoly = hPoly + Hk(i,[0 1 2]);
end
hPoly

hPoly = 1×3

-1 1 5

Matlab hermiteH documentation.

The code is Not right.

1. The MATLAB hermiteH is the "physicist's Hermite polynomials" (i.e. ) and not the "probabilist's Hermite polynomials" (i.e. ) I'm looking for but I "fixed" it with the extra "" term, but I'm not sure.
2. The result is wrong and I cant get it right.
3. The end goal is to implement the Wolfram syntax into a functional MATLAB code with k and x variables. note: to make multiple polynomials in one vector with the hermiteH command it requires k to be a vector (i.e. 1:7) and x to be symbolic so basically k is a vector and x is a vector and make the implementation versetile.
4. The problem I think is the first term in the result adds the first Hermite k times, the second Hermite k-1 times and so on...

Any advice on the matter is much obliged.

[Dyuman Joshi]

I have derived this from the code I wrote for "physicist's Hermite polynomials" and transformed it according to the rescaling equations from Wikipedia.

The code gives the coefficients of polynomials in the form of a vector. However, this code is susceptible to rounding errors as it includes multiplication with irrational values.

I'll update the code with the actual algorithm for probabilist's Hermite polynomial when I'll get the time. Hope this helps.

n=5;
%physicist's hermite polynomial
q=1;
for i=1:n
    y=conv([-2 0],q);
    der=polyder(q);
    q=y+[zeros(1,numel(y)-numel(der)) der];
end
y=sqrt(2);
q=q*(-1)^n

q = 1×6

32 0 -160 0 120 0

%probabilit's hermite polynomial
p=q.*y^(-n).*y.^-(n:-1:0)

p = 1×6

1.0000 0 -10.0000 0 15.0000 0

[Alex Rudyak]

1. I didn't quite understand where did the in the conv function came from.
2. Is this a one to one implementation of the equation?
3. I'm not familiar with polyder I'll look into it

Thank you for this, I'll have to try it out first

[John D'Errico]

@Alex Rudyak

conv is the classic MATLAB trick for multiplication of polynomials or many other things. You can use it to build polynomials, etc. I even used it in my VPI toolbox to perform multiplication of variable precision numbers. For example, suppose you want to multiply a pair of polynomials? I'll do it using syms first.

Suppose I want to multiply (x+1) by (x^2 - 3*x - 2)? Firt, in symbolic form so we can compare the results.

syms x
expand((x+1)*(x^2 - 3*x - 2))

ans = 

But now recognize that these polynomials can be represented as a list of coefficients in standard MATLAB form.

conv([1,1],[1 -3 -2])

ans = 1×4

1 -2 -5 -2

Do you see the multiplication worked perfectly? As a result, we now have a list of coefficients of the product polynomial. Try it with other products of polynomials, and it always works, at least as long as there are no problems with the precision of the coefficients in double precision arithmetic.

Why does it work? Think about what conv does. Expand it out, and you will see that convolution does exactly what you want with the product of those polynomials. (Conv is a tremendously valuable tool, useful in many areas of mathematics. I cannot begin to count the number of ways I have used it over the years.)

[Alex Rudyak]

@John D'Errico I've completly missed this one, thank you it's clearer now, I'll note it down for future use

Answer 1

[John D'Errico]

First, the two sets of polynomials are effectively the same thing, but for a transformation. As you said, the two are just orthogonal with respect to a subtly different version of the same exponential form. Here, for example, is the degree 3 polynomial, and that matches what wikipedia shows. First, note the normalization. These are not normalized to have the coefficient of the highest degree term 1.

syms x
P3 = hermiteH(3,x)

P3 = 

This is what Wikipedia refers to as a physicist's polynomials. They are orthognal wrt exp(-x^2). So effectively derived from an erf, which makes sense in context. For example...

P2 = hermiteH(2,x)

P2 = 

then we would expect to find this integral to be zero.

int(P2*P3*exp(-x^2),[-inf,inf])

ans = 

0

Every once in a while something works. ;-) But not everything. I prefer my orthognal polynomials to be normalized to have unit norm too. And that would mean this next would yield 1. Such is life. Just a scale factor, and I can live with it.

int(P3*P3*exp(-x^2),[-inf,inf])

ans = 

The Probabilists form will be orthogonal to exp(-x^2/2). We should recognize this from the Gaussian (Normal) distribution PDF. And we can get one from the other.

S2 = sqrt(sym(2));
hermiteH(3,x/S2)

ans = 

But we need to normalize it, so the highest order term has coefficient 1, at least if you want to match the table. And to do that, you just divide by 2^(n/2), So I'll write a little function that does it:

hermiteProb = @(x,n) simplify(hermiteH(n,x/S2)/S2^n);

P2prob = hermiteProb(x,2)

P2prob = 

P3prob = hermiteProb(x,3)

P3prob = 

int(P2prob*P3prob*exp(-x^2/2),[-inf,inf])

ans = 

0

Again, they are orthogonal wrt the indicated weight function. (WHEW! Got lucky again!) Though they are not normalized so the integral is 1, as would be expected if they were orthonormal. Oh well, I guess nothing is ever perfect. Still they have the form you desire.

int(P3*P3*exp(-x^2/2),[-inf,inf])

ans = 

Still not my favorite normalization. Oh well.

Finally, if you want to evaluate these polynomials at a list of points, that is pretty easy. You could convert them to a function, using matlabFunction.

P3probfun = matlabFunction(P3prob)

P3probfun = function_handle with value:
    @(x)x.*-3.0+x.^3

P3probfun(-1.5:.5:1.5)

ans = 1×7

1.1250 2.0000 1.3750 0 -1.3750 -2.0000 -1.1250

Or plot it.

X = linspace(-5,5);
plot(X,P3probfun(X),'-')
grid on

[Alex Rudyak]

Simplify is presenting an error, I've removed it and S2 term and it works:

Hk = @(x,n) hermiteH(n,x/sqrt(2))/sqrt(2)^n;

Hk([0 1 2 4 20],2)

ans = 1×5

-1.0000 0 3.0000 15.0000 399.0000

I'll try and use the

% P3probfun = matlabFunction(P3prob)

example you provided to use a for loop.

As thorough as this explanation and helpful, I'm still left with a need to run a vector through 7 different polynomials and add them up, I'm currently working on an implementation based on this answer

Thank you for your time, I'll update on the results

update:

Perhaps something like this:

Hk = @(x,n) hermiteH(n,x/sqrt(2))/sqrt(2)^n;
x = [-0.5 -0.1 0 0.1 0.5 1 2 3];
Hkvec = Hk(x,1) + Hk(x,2) + Hk(x,3) + Hk(x,4) + Hk(x,5) + Hk(x,6) + Hk(x,7)

Hkvec = 1×8

30.7578 -3.4972 -13.0000 -21.7056 -38.4766 -1.0000 59.0000 -415.0000

Since I wont be taking more then the 7th this might be valid, I'll keep updating with news

[John D'Errico]

You need to define S2 in advance for that part to work. If you do, then when you create the function handle, S2 gets encapsulated into the function handle workspace, and it works.

As for writing a loop to evaluate 7 different polynomials and sum them in some way, that is your job. You should know how to use a loop. Its just a loop after all. You can store the various polynomials in a cell array if you want, then evaluate them. If you store the coefficients, then you can use polyval to evaluate them.

Regardless, I think you now know what you need to know, and that is what matters.

[Alex Rudyak]

Thank you, much obliged!

[Alex Rudyak]

Hey just a quick update,

I've calculated the polynomials and just matlab functioned it all like so:

syms x
S2 = sqrt(sym(2));
hermiteProb = @(x,n) simplify(hermiteH(n,x/S2)/S2^n);

P0prob = hermiteProb(x,0);
P1prob = hermiteProb(x,1);
P2prob = hermiteProb(x,2);
P3prob = hermiteProb(x,3);
P4prob = hermiteProb(x,4);
P5prob = hermiteProb(x,5);
P6prob = hermiteProb(x,6);
P7prob = hermiteProb(x,7);
PNprobVec = matlabFunction([P0prob P1prob P2prob P3prob P4prob P5prob P6prob P7prob]);
x = [-0.5 -0.1 0 0.1 0.5 1 2 3];
g = PNprobVec(x);
g = reshape(g(2:end),8,7);
PNprobMat = table(g(:,1), g(:,2), g(:,3), g(:,4), g(:,5), g(:,6), g(:,7));
PNprobMat.Properties.VariableNames = ["P1 Answer","P2 Answer","P3 Answer","P4 Answer","P5 Answer","P6 Answer","P7 Answer"]

PNprobMat = 8×7 table
    P1 Answer    P2 Answer    P3 Answer    P4 Answer    P5 Answer    P6 Answer    P7 Answer
    _________    _________    _________    _________    _________    _________    _________

      -0.5         -0.75        1.375       1.5625       -6.2812      -4.6719       40.023 
      -0.1         -0.99        0.299       2.9401         -1.49      -14.551       10.395 
         0            -1            0            3             0          -15            0 
       0.1         -0.99       -0.299       2.9401          1.49      -14.551      -10.395 
       0.5         -0.75       -1.375       1.5625        6.2812      -4.6719      -40.023 
         1             0           -2           -2             6           16          -20 
         2             3            2           -5           -18          -11           86 
         3             8           18           30            18          -96         -396 

Now I need to apply an integration on the matrix, but thats a question for another day.

Thank you for your time!