Hermite polynomial method, incorrect coefficients

I want to obtain coefficients through Hermite polynomial method but the coefficients that I obtained trough my code are incorrect

clear variables;
close all;
clc;

x = [0 0.5];
y = [1 2.71828]; 
yp = [2 5.43656]; 

n = length(x)-1;
m = 2; 

x_bar = zeros((n+1)*(m+1),1);
H = zeros((n+1)*(m+1),(n+1)*(m+1));

for k = 0:n
    j1 = 1:m+1;
    j2 = 2:m+1;


    x_bar((m+1)*k+j1) = x(k+1);
    H((m+1)*k+j1,1) = y(k+1);
    H((m+1)*k+j2,2) = yp(k+1);


    if k~=0
    H((m+1)*k+1,2) = (H((m+1)*k+1,1)-H((m+1)*k,1))/(x_bar((m+1)*k+1)-x_bar((m+1)*k));
    end
end

idx = [4 5];

for i=1:length(idx)
    for j = 3
        H(idx(i),j) = (H(idx(i),j-1)-H((idx(i)-1),j-1))/(x_bar(idx(i))-x_bar(idx(i)-j+1));
    end
end


for i = 4:(n+1)*(m+1)
    for j = 4:i
        H(i,j) = (H(i,j-1)-H(i-1,j-1))/(x_bar(i)-x_bar(i-j+1));
    end
end


x_bar;
H

coeff = diag(H)

[Dyuman Joshi]

Which type of Hermite polynomial are you calculating? Probablist's or Physicist's?

[David Cuenca]

Probabilistic mate

Answer 1

[John D'Errico]

I could swear I answered this question only recently. Maybe not for you, and I don't feel like a search through the many thousands of my answers. But you have made at least a serious effort, and your code loooks like it is at least going in the right direction.

A quick look at your code shows some serious problems though. First, don't code the value exp(1) as a number. It already exists in MATLAB.

format long g 
exp(1)

ans = 
          2.71828182845905

And, honestly, trying to get through the completely uncommented code you have written will be almost impossible. Sorry. That is what happens when you write code like that. Others don't want to read it.

So what order polynomial are you trying to generate? It looks like you want to generate a complete table of polynomials, so this is surely homework. And the fact that you defined exp(1) as 2.71828 is another clue.

The probabilist form is different only from the alternate physicist form in that they are orthogonal with respect to different weight functions. So the physicist form (H_n) is orthogonal to

exp(-x^2)

while the probabilist form polynomials (He_n) are orthogonal to

exp(-x^2/2)/sqrt(2*pi)

effectively this results in nothing more than a simple scale transformation of the polynomials. But either is easy to generate. Use a recurrence relation. For the probabilist ones, that would be:

He_(n+1) = x*He_n - diff(He_n)

In both cases, the zero-th order polynomial is the constant 1. So first, using symbolic tools, I might do this:

M = 5; % highest order
syms x
He = zeros(M+1,1,'sym'); % offset by 1, because MATLAB has index origin 1
He(1) = sym(1);
for m = 1:M
    He(m+1) = expand(x*He(m) - diff(He(m),x,1)); % recurrence relation
end
He

He = 

These are of course the same as those seen in the table on the wiki page.

We can even verify they are orthogonal, wrt exp(-x^2/2). This next should be zero.

int(He(2)*He(3)*exp(-x^2/2)/sqrt(2*sym(pi)),[-inf,inf])

ans = 

0

Are they orthonormal? If so, this next should be 1.

int(He(2)*He(2)*exp(-x^2/2)/sqrt(2*sym(pi)),[-inf,inf])

ans = 

1

So they have been normalized where the dot product is 1. So these polynomials are orthonormal wrt that dot product. (This is not a general proof of course.)

Doing the same thing in an array, where we would generate a list of coefficients is trivial of course. Just use the same recurrence relation.

Hecoeff = zeros(m+1,m+1);
Hecoeff(1,1) = 1;
for m = 1:M
  Hecoeff(m+1,1:m+1) = [0,Hecoeff(m,1:m)] - [(1:m-1).*Hecoeff(m,2:m),0 0]; % recurrence relation
end
Hecoeff

Hecoeff = 6×6

1 0 0 0 0 0 0 1 0 0 0 0 -1 0 1 0 0 0 0 -3 0 1 0 0 3 0 -6 0 1 0 0 15 0 -10 0 1

If you inspect these polynomial coefficients, you will see they are ordered so the lowest order term is first. This is of course backwards form the way polyval would expect. So, sue me. But they are in fact the same polynomials as contained in the array He. We can even recover them in symbolic form as:

Hecoeff*x.^(0:M).'

ans = 

Of course, we can evaluate those polynomials easily enough. For example, we can plot them.

fplot(He(4)) % the cubic polynomial

Or I might evaluate the 4'th degree polynomial at some point. If you use polyval, don't forget to flip the coefficients.

polyval(flip(Hecoeff(5,1:5)),1.75)

ans = 
               -5.99609375