How can I normalize a function solved numerically?

I'm trying to make the pdf of a quantum harmonic oscillator, once the equation is solved numerically, the pdf area should be 1 but instead it is 2.2604e-28. Thanks for advanced

function test(n,y0,dy0)
syms y(x)
m=9.1093821545*1e-31; k=m*(4.57*1e14)^2;  w=sqrt(k/m); hb=(6.6260695729*1e-34)/(2*pi);  En=hb*w*(n+.5); c1=-2*m/(hb^2); c_i0=sqrt(hb/(m*w));
[V] = odeToVectorField(diff(y,x,2) == c1*(En-.5*k*x.^2).*y);
x1=0;
x2=2e-9;
M = matlabFunction(V,'vars',{'x','Y'});
[p,q]= ode45(M,[x1 x2],[y0 dy0]);
dens=(q(:,1).*q(:,1));
% A plotear
plot([-p p],En+dens,'-k') %psi quadrat
%plot(-p,En+q(:,1),'-k') psi
hold on
%plot(p,En+q(:,1),'-k') psi
potencial=0.5*k*p.^2;
plot([-p p],potencial,'-r') %potencial
plot([-p p], ones(1,length(p))*En,'-b') %En
% natural frequency of the oscillator w = 4.57e14 Hz
hold off
grid on
u=zeros(1,length(dens));
2*trapz(p,dens) %area of the pdf
end

[Star Strider]

What are typical values for ‘n’, ‘y0’, and ‘dy0’ ?

[david feldstein]

they're initial conditions, i.e: 0,5.0331e-10,0

[Walter Roberson]

Instead of doing the trapz for dens over time, you could add another variable to the ode which is V(1).^2 . The output at the corresponding position would automatically be the integral of that value.

[david feldstein]

thanks! but do you know if there's some method to normalize the function numerically? or it can only be done analytically?

[Torsten]

Walter means that you should add the equation

dy(3)/dt = y(1)^2

to your system of ordinary differential equations, integrate numerically and afterwards normalize y(1) according to

dense=dense/q(end,3).

Best wishes

Torsten.

[david feldstein]

right thanks!

Answer 1

[David Goodmanson]

Hi david,

A normalization integral that small shows that the scaling could be improved upon, but as a workaround is this what you are looking for?

densnorm = dens/(2*trapz(p,dens));
2*trapz(p,densnorm)                     %area of the pdf
ans = 1

[david feldstein]

thanks David! the plot would be this one right?

plot([-p p],densnorm)

[David Goodmanson]

yes, that's the new one to use. Maybe I should have just said dens = dens / (2*trapz(p,dens)) but I wanted to emphasize normalization.

Incidentally, if you go to rescaled units for x as you almost did, then things are easier:

xs = sqrt(hb/(m*w))   length scale, same as your c_i0
Es = hb*w             energy scale
Fn = En/Es = n+1/2    normalized energy
u = x/xs              normalized x

In that case all the constants disappear from the differential equation, leaving just the factors of 1/2 and you can solve

[V] = odeToVectorField( (-1/2)*diff(y,u,2) == (Fn-(1/2)*u.^2).*y );

The integration limits on u are something like 0 to 10.

To normalize the resulting density with the new p array (which represents u) you do exactly like before, dens = dens / (2*trapz(p,dens)). It's easier to do plots since now the potential energy is (1/2)*p.^2 and fits on the plot without rescaling.

Then if you want to get the real distances back again, since p means u right now, you would use x = p*xs as the new variable. If you wanted a normalized density in terms of x you would invoke (guess what) dens = dens / (2*trapz(x,dens)).

[david feldstein]

thank you! I'm definitely trying your method