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
6 个评论
Star Strider
2017-11-30
What are typical values for ‘n’, ‘y0’, and ‘dy0’ ?
david feldstein
2017-11-30
Walter Roberson
2017-12-1
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
2017-12-1
Torsten
2017-12-1
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
2017-12-1
采纳的回答
David Goodmanson
2017-12-1
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
3 个评论
david feldstein
2017-12-1
David Goodmanson
2017-12-2
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
2017-12-3
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Quantum Mechanics 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
