Please give us the details on what you are trying to integrate, and over what range.
Numerical Integration in matlab
2 次查看（过去 30 天）
显示 更早的评论
what is reliable way to perform a integartion numerically in matlab? I would like to get numerical result after integration of a large function .
Is that Trapezoidal rule or Simpson's rule is reliable to do that? Pl somebody tell me what is the best way for that.
Otherwise, when I use 'int()' command directly, Matlab takes very very huge time almost 5-6 hours.
4 个评论
AVM
2020-1-22
@James: Thanks . Pl see this.
clc
close all
syms 'theta' 'phi' g
g=1;
thet=4;
h=[ cos(theta) 0 sin(theta)*exp(-1i*phi) 0
0 cos(theta) g sin(theta)*exp(-1i*phi)
sin(theta)*exp(1i*phi) g -cos(theta) 0
0 sin(theta)*exp(1i*phi) 0 -cos(theta) ];
a1=(exp(-phi*2i)*(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2)^(3/2))/(g*sin(theta)^2) - (exp(-phi*2i)*cos(theta)*(g^2 + cos(theta)^2 + sin(theta)^2))/(g*sin(theta)^2) + (exp(-phi*2i)*cos(theta)*(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2))/(g*sin(theta)^2) - (exp(-phi*2i)*(g^2 + cos(theta)^2 + sin(theta)^2)*(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2)^(1/2))/(g*sin(theta)^2);
b1=(exp(-phi*1i)*(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2)^(1/2))/sin(theta) + (exp(-phi*1i)*cos(theta))/sin(theta);
c1=-(exp(-phi*1i)*(cos(theta)^2 + sin(theta)^2))/(g*sin(theta)) + (exp(-phi*1i)*(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2))/(g*sin(theta));
d1=1;
m=sqrt(a1*conj(a1)+b1*conj(b1)+c1*conj(c1)+d1*conj(d1));
u=1/m*[ a1
b1
c1
d1 ];
v=diff(u,phi);
r=dot(u,v);
f(theta,g)=1i/pi*int(r,phi,0,2*pi); %% This integration taking so much time
f = simplify(f, 'Steps',500);
ffcn = matlabFunction(f);
theta = linspace(0.001,4, 30);
g = linspace(0.001,10, 30);
[Th,G] = meshgrid(theta, g);
F=ffcn(Th,G);
%max(imag(F(:)))
%min(imag(F(:)))
figure
mesh(Th, G, F)
colormap(cool)
grid on
xlabel('\bf\theta','FontSize',14)
ylabel('\bf\alpha','FontSize',14)
zlabel('\bf\itf','FontSize',14)
%% This entire process takes long time, I quit that in between.
Walter Roberson
2020-1-22
编辑：Walter Roberson
2020-1-22
simplify( r) before doing the integration; it compacts down quite a bit. Also,
assume(phi>=0)
assume(theta>=0)
before doing the simplify()
采纳的回答
Walter Roberson
2020-1-22
There is no reliable numeric integration method in any programming language. For any given numeric integration method, it is possible to construct a function which the numeric integration method will return an arbitrarily wrong solution.
Trapazoidal Rule and Simpson's Rule are not reliable at all.
You could consider using integral(); provide AbsTol and RelTol parameters, and provide WayPoints whenever meaningful.
You could consider using vpaintegral() with similar parameters.
23 个评论
Walter Roberson
2020-1-22
Once you
assume(phi >= 0)
which we can do because you integrate r over phi from 0 to 2*Pi
then when you simplify(r ) the phi drops out and you end up with an expression that is only in theta. Because it is constant in the variable of integration, the integral is just the difference in integration bounds times the value of the simplified expression.
AVM
2020-1-22
@ Walter: Can I take steps size smaller than 500? Bcz it is taking so much time for simplification r.
Walter Roberson
2020-1-22
编辑：Walter Roberson
2020-1-22
>> assume(phi>=0); assume(theta>=0);
>> tic;symvar(simplify(r)),toc
ans =
theta
Elapsed time is 0.389983 seconds.
And therefore there is no need to int() with respect to phi, since phi dropped out of the expression.
AVM
2020-1-22
@Walter: Thanks a lot. The plotting is beautifully coming. I am really greatful to you.
AVM
2020-1-22
@Walter: Is there any way to call the particular eigen vector of a matrix? Actually, in my current programming, this a1,b1,c1,d1 are basically components of a particular eigen vector of h matirx. I just paste them from command palte. But I would like to avoid that big expression. Is there any way to call it directly?
Walter Roberson
2020-1-22
Eigenvectors are not unique so it is not clear what a particular eigenvector is.
If you use the two output form of eig() then some of the eigenvectors appear as columns of the first output. If the eigenvalues are unique then the eigenvectors output span the space.
AVM
2020-1-23
编辑：AVM
2020-1-23
@walter: Thanks ..that I know. Here I attach my code, Why this is taking so much time to reveal the eigen values of the finalrho? As if the process is never ending. You Pl look my code.
clc
close all
syms 'theta' 'alpha' 'phi' k
%% sigma matrices are here
sigma1=[0 1;1 0];
sigma2=[0 -1i;1i 0];
sigma3=[1 0;0 -1];
sigmap=1/2*(sigma1+1i*sigma2);
sigmam=1/2*(sigma1-1i*sigma2);
%% unit vetor on sphere
n=[sin(theta)*cos(phi);sin(theta)*sin(phi);cos(theta)];
identity1=eye(2);
identity2=identity1;
p1=sigma1*sin(theta)*cos(phi)+sigma2*sin(theta)*sin(phi)+sigma3*cos(theta);
p2=kron(sigmap,sigmap);
p3=kron(sigmam,sigmam);
h=1/2*alpha*kron(p1,identity2)+k*(p2+p3);
[V,L]=eig(h);
%% The componens of one of the eigen vectors of h above pasted from command palte
a1=(4*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4)^(3/2))/(alpha^2*k*cos(phi)^2*sin(theta)^2 + alpha^2*k*sin(phi)^2*sin(theta)^2) - (4*k^2*cos(theta) + alpha^2*cos(theta)^3 + alpha^2*cos(phi)^2*cos(theta)*sin(theta)^2 + alpha^2*cos(theta)*sin(phi)^2*sin(theta)^2)/(2*(alpha*k*cos(phi)^2*sin(theta)^2 + alpha*k*sin(phi)^2*sin(theta)^2)) + (2*cos(theta)*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4))/(alpha*k*cos(phi)^2*sin(theta)^2 + alpha*k*sin(phi)^2*sin(theta)^2) - ((alpha^2*cos(theta)^2 + 4*k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4)^(1/2))/(alpha^2*k*cos(phi)^2*sin(theta)^2 + alpha^2*k*sin(phi)^2*sin(theta)^2);
b1=-(((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4)^(3/2)*8i)/(alpha^3*cos(phi)^3*sin(theta)^3*1i - alpha^3*sin(phi)^3*sin(theta)^3 + alpha^3*cos(phi)*sin(phi)^2*sin(theta)^3*1i - alpha^3*cos(phi)^2*sin(phi)*sin(theta)^3) + (cos(theta)*(alpha^2*cos(theta)^2 + 4*k^2 + 2*alpha^2*cos(phi)^2*sin(theta)^2 + 2*alpha^2*sin(phi)^2*sin(theta)^2))/(alpha^2*cos(phi)^3*sin(theta)^3 + alpha^2*sin(phi)^3*sin(theta)^3*1i + alpha^2*cos(phi)*sin(phi)^2*sin(theta)^3 + alpha^2*cos(phi)^2*sin(phi)*sin(theta)^3*1i) + (2*(alpha^2*cos(theta)^2 + 4*k^2 + 2*alpha^2*cos(phi)^2*sin(theta)^2 + 2*alpha^2*sin(phi)^2*sin(theta)^2)*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4)^(1/2))/(alpha^3*cos(phi)^3*sin(theta)^3 + alpha^3*sin(phi)^3*sin(theta)^3*1i + alpha^3*cos(phi)*sin(phi)^2*sin(theta)^3 + alpha^3*cos(phi)^2*sin(phi)*sin(theta)^3*1i) - (cos(theta)*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4)*4i)/(alpha^2*cos(phi)^3*sin(theta)^3*1i - alpha^2*sin(phi)^3*sin(theta)^3 + alpha^2*cos(phi)*sin(phi)^2*sin(theta)^3*1i - alpha^2*cos(phi)^2*sin(phi)*sin(theta)^3);
c1=-((alpha^2*cos(theta)^2)/2 + 2*k^2 + (alpha^2*cos(phi)^2*sin(theta)^2)/2 + (alpha^2*sin(phi)^2*sin(theta)^2)/2)/(2*((alpha*k*cos(phi)*sin(theta))/2 - (alpha*k*sin(phi)*sin(theta)*1i)/2)) + (2*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4))/(alpha*k*cos(phi)*sin(theta) - alpha*k*sin(phi)*sin(theta)*1i);
d1=1;
m=sqrt(a1*conj(a1)+b1*conj(b1)+c1*conj(c1)+d1*conj(d1));
u=1/m*[a1;b1;c1;d1];
rho=u*ctranspose(u);
y=kron(sigma2,sigma2);
rhofinal=rho*y*ctranspose(rho)*y;
q=eig(rhofinal)
Moreover, is there any way to call this eigen vector for next operation without expressing its large form. Here a1,b1,c1 and d1 are basically components of a particular eigen vectors of matrix h. I pasted it from command plate
Walter Roberson
2020-1-23
If I respond, are you going to delete my contributions later, like you did in the 3d plotting question?
AVM
2020-1-23
@walter: Thanks for replying me.
Apologize for my nonsense activities.. I am really new user this matlab. I thought that as I asked you people so many times regarding my corrections and rectifications, so better to clear my all these garbage kind of thing so that some important and summary conversation remain. But I really don't know that it actually causes some harmful to you.
Pl help me first to recover all that deleted comment by me. I am really feeling a guilty.
AVM
2020-1-23
@walter: PL help me to recover that conversation first if possible. Otherwise, I always feel ashame to ask you further.
AVM
2020-1-23
@walter : I am really sorry. This happens beyond my knowldege. I am promising you that next time such kind of nonsense would not happen.
I must be carefull about all these options. Pl pardon me if possible.
AVM
2020-1-24
编辑：AVM
2020-1-24
@walter : Pl help me to solve the issue.
I would like to get eigenvalues of matrix in my code given below, the matlab taking so much time as if never ending processws. Pl somebody hepl me how to find the eigenvalues quickly.
clc
close all
syms 'theta' 'alpha' 'phi' k
%% sigma matrices are here
sigma1=[0 1;1 0];
sigma2=[0 -1i;1i 0];
sigma3=[1 0;0 -1];
sigmap=1/2*(sigma1+1i*sigma2);
sigmam=1/2*(sigma1-1i*sigma2);
%% unit vetor on sphere
n=[sin(theta)*cos(phi);sin(theta)*sin(phi);cos(theta)];
identity1=eye(2);
identity2=identity1;
p1=sigma1*sin(theta)*cos(phi)+sigma2*sin(theta)*sin(phi)+sigma3*cos(theta);
p2=kron(sigmap,sigmap);
p3=kron(sigmam,sigmam);
h=1/2*alpha*kron(p1,identity2)+k*(p2+p3);
[V,L]=eig(h);
%% The componens of one of the eigen vectors of h above pasted from command palte
a1=(4*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4)^(3/2))/(alpha^2*k*cos(phi)^2*sin(theta)^2 + alpha^2*k*sin(phi)^2*sin(theta)^2) - (4*k^2*cos(theta) + alpha^2*cos(theta)^3 + alpha^2*cos(phi)^2*cos(theta)*sin(theta)^2 + alpha^2*cos(theta)*sin(phi)^2*sin(theta)^2)/(2*(alpha*k*cos(phi)^2*sin(theta)^2 + alpha*k*sin(phi)^2*sin(theta)^2)) + (2*cos(theta)*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4))/(alpha*k*cos(phi)^2*sin(theta)^2 + alpha*k*sin(phi)^2*sin(theta)^2) - ((alpha^2*cos(theta)^2 + 4*k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4)^(1/2))/(alpha^2*k*cos(phi)^2*sin(theta)^2 + alpha^2*k*sin(phi)^2*sin(theta)^2);
b1=-(((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4)^(3/2)*8i)/(alpha^3*cos(phi)^3*sin(theta)^3*1i - alpha^3*sin(phi)^3*sin(theta)^3 + alpha^3*cos(phi)*sin(phi)^2*sin(theta)^3*1i - alpha^3*cos(phi)^2*sin(phi)*sin(theta)^3) + (cos(theta)*(alpha^2*cos(theta)^2 + 4*k^2 + 2*alpha^2*cos(phi)^2*sin(theta)^2 + 2*alpha^2*sin(phi)^2*sin(theta)^2))/(alpha^2*cos(phi)^3*sin(theta)^3 + alpha^2*sin(phi)^3*sin(theta)^3*1i + alpha^2*cos(phi)*sin(phi)^2*sin(theta)^3 + alpha^2*cos(phi)^2*sin(phi)*sin(theta)^3*1i) + (2*(alpha^2*cos(theta)^2 + 4*k^2 + 2*alpha^2*cos(phi)^2*sin(theta)^2 + 2*alpha^2*sin(phi)^2*sin(theta)^2)*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4)^(1/2))/(alpha^3*cos(phi)^3*sin(theta)^3 + alpha^3*sin(phi)^3*sin(theta)^3*1i + alpha^3*cos(phi)*sin(phi)^2*sin(theta)^3 + alpha^3*cos(phi)^2*sin(phi)*sin(theta)^3*1i) - (cos(theta)*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4)*4i)/(alpha^2*cos(phi)^3*sin(theta)^3*1i - alpha^2*sin(phi)^3*sin(theta)^3 + alpha^2*cos(phi)*sin(phi)^2*sin(theta)^3*1i - alpha^2*cos(phi)^2*sin(phi)*sin(theta)^3);
c1=-((alpha^2*cos(theta)^2)/2 + 2*k^2 + (alpha^2*cos(phi)^2*sin(theta)^2)/2 + (alpha^2*sin(phi)^2*sin(theta)^2)/2)/(2*((alpha*k*cos(phi)*sin(theta))/2 - (alpha*k*sin(phi)*sin(theta)*1i)/2)) + (2*((alpha^2*cos(theta)^2)/4 - (k*(k^2 + alpha^2*cos(phi)^2*sin(theta)^2 + alpha^2*sin(phi)^2*sin(theta)^2)^(1/2))/2 + k^2/2 + (alpha^2*cos(phi)^2*sin(theta)^2)/4 + (alpha^2*sin(phi)^2*sin(theta)^2)/4))/(alpha*k*cos(phi)*sin(theta) - alpha*k*sin(phi)*sin(theta)*1i);
d1=1;
m=sqrt(a1*conj(a1)+b1*conj(b1)+c1*conj(c1)+d1*conj(d1));
u=1/m*[a1;b1;c1;d1];
rho=u*ctranspose(u);
y=kron(sigma2,sigma2);
rhofinal=rho*y*ctranspose(rho)*y;
q=eig(rhofinal)
Second thing , is there any way to call the partulcar eigen vector of a given matirx. Because, here by the expression of a1,b1,c1 and d1 I actually have written the components of one of the eigenvector of h from command plate. But I think it is not an efficient way to write.
Pl help to do that.
Walter Roberson
2020-1-24
"is there any way to call the partulcar eigen vector of a given matirx."
That question is not meaningful as phrased.
Use the two output form of eig() and the eigenvectors will be the columns of the first output. You can index them directly by column number.
AVM
2020-1-24
@Walter: Here I try to calculate the concurrenec for qubit system.
Concurrence ©= max{0, e1-e2-e3-e4} where e1,e2,e3,e4 are the eigenvalues of the density matrix of R, the form given in my code.
Here e1>=e2>=e3>=e4.
Pl help me to get the expression for concurrence.
clc
close all
syms 'theta' 'alpha' 'phi' k
%% sigma matrices are here
sigma1=[0 1;1 0];
sigma2=[0 -1i;1i 0];
sigma3=[1 0;0 -1];
sigmap=1/2*(sigma1+1i*sigma2);
sigmam=1/2*(sigma1-1i*sigma2);
%% unit vetor on sphere
n=[sin(theta)*cos(phi);sin(theta)*sin(phi);cos(theta)];
identity1=eye(2);
identity2=identity1;
p1=sigma1*sin(theta)*cos(phi)+sigma2*sin(theta)*sin(phi)+sigma3*cos(theta);
p2=kron(sigmap,sigmap);
p3=kron(sigmam,sigmam);
h=1/2*alpha*kron(p1,identity2)+k*(p2+p3);
[V,L]=eig(h);
%%state vector
x1=V(:,1);
m=sqrt(dot(x1,x1));
%% Normalized state vector
x2=x1/m;
%% Density matrix
rho=x2*ctranspose(x2);
y=kron(sigma2,sigma2);
function c=concurrence(rho);
R=rho*y*ctranspose(rho)*y;
% Real is needed since MATLAB always adds a small imaginary part ...
e=real(sqrt(eig(R)));
e=-sort(-e);
c=max(e(1)-e(2)-e(3)-e(4),0)
end
I need the expression of c, but here nothing is coming out.
AVM
2020-1-25
@walter: How to resume the Matlab programming when it is pause situation?? And how to stop it when it is showing busy below without closing Matlab software itself?
Walter Roberson
2020-1-26
What is a "pause situation" ?
If you are stopped in the debugger and wish to continue then
dbcont
If it is busy and you have the editor open, in more recent versions you can click on the Pause button . It might take some time before you get the command line control back.
If the delay for Pause is getting too long, you can control-C in the command window. There are some circumstances under which it can take a fair time to get the command line control back anyhow.
If you managed to use up all system memory and you are swapping to disk, then most often you should use your operating system task killer to kill MATLAB .
AVM
2020-1-27
Thanks...Are the latest version of matlab little bit faster than 2018a? I am using 2018a. It is taking so much time for even 'simplify(f)' kind of operation and also for any complicated mathematical work.
Walter Roberson
2020-1-27
编辑：Walter Roberson
2020-1-27
There is no faster version of symbolic simplification, except possibly some fairly small gains in performance in newer versions.
Symbolic work is done using compiled libraries that are not written in MATLAB itself, so improvements in the MATLAB Execution Engine that have been made do not improve the symbolic engine.
If you are doing numeric integration on a system that has only one free variable, then use vpaintegral() instead of int() -- but integral() will likely be faster than vpaintegral()
更多回答（2 个）
AVM
2020-1-26
@walter: How do I can normalize an eigen vector in Matlab? PL tell me what is the corresponding command for the normalization?
4 个评论
Walter Roberson
2020-1-27
编辑：Walter Roberson
2020-1-27
V_normalized = V ./ sqrt(sum(V.^2)); %requires R2016b or newer
AVM
2020-1-27
@walter: With this normalized command when I try to run the following progamming, it is taking so much time as near to 1 hour and the pc start to hang frequently. I had forcefully quit that. You, pl have a look on my code and suugest me what to do. I am using matlab R2018a version.
clc;clear;
syms theta phi g
%% h is 4*4 matrix
h=[ cos(theta) 0 sin(theta)*exp(-1i*phi) g
0 cos(theta) 0 sin(theta)*exp(-1i*phi)
sin(theta)*exp(1i*phi) 0 -cos(theta) 0
g sin(theta)*exp(1i*phi) 0 -cos(theta) ];
%% a1,b1,c1 and d1 are the components of first eigenvector of h; these are pasted from command plate
%a1=(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2)^(3/2)/(g*sin(theta)^2) - (cos(theta)^3 + cos(theta)*sin(theta)^2 + g^2*cos(theta))/(g*sin(theta)^2) + (cos(theta)*(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2))/(g*sin(theta)^2) - ((g^2 + cos(theta)^2 + sin(theta)^2)*(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2)^(1/2))/(g*sin(theta)^2);
%b1=-(exp(-phi*1i)*(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2)^(3/2))/sin(theta)^3 + (exp(-phi*1i)*(cos(theta)^2 + 2*sin(theta)^2 + g^2)*(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2)^(1/2))/sin(theta)^3 - (exp(-phi*1i)*cos(theta)*(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2))/sin(theta)^3 + (exp(-phi*1i)*cos(theta)*(cos(theta)^2 + 2*sin(theta)^2 + g^2))/sin(theta)^3;
%c1= -(g^2*exp(phi*1i) + exp(phi*1i)*cos(theta)^2 + exp(phi*1i)*sin(theta)^2)/(g*sin(theta)) + (exp(phi*1i)*(cos(theta)^2 + sin(theta)^2 + g^2/2 - (g*(4*sin(theta)^2 + g^2)^(1/2))/2))/(g*sin(theta));
%d1=1;
%m=sqrt(a1*conj(a1)+b1*conj(b1)+c1*conj(c1)+d1*conj(d1));
%% normalized first eigen vector of h
%u=1/m*[a1;b1;c1;d1];
%% differently
[V,L]=eig(h);
u=V(:,1)./sqrt(sum(V(:,1).^2)); %% this call aslo for normalised eigenvector of h
x=diff(u,phi);
r=dot(u,x);
assume(theta>=0);
assume(phi>=0);
r=simplify(r,'Steps',100);
f(theta,g)=1i/pi*int(r,phi,0,2*pi);
f=simplify(f,'Steps',100);
ffcn = matlabFunction(f);
theta = linspace(0.001,4, 30);
g = linspace(.001,2, 30);
[Th,G] = ndgrid(theta, g);
%F=abs(ffcn(Th,Al));
F=ffcn(Th,G);
%max(imag(F(:)))
%min(imag(F(:)))
figure
meshc(Th,G, F)
colormap(cool)
grid on
xlabel('\bf\theta','FontSize',14)
ylabel('\bf\alpha','FontSize',14)
zlabel('\bf\itf','FontSize',14)
But when I normalize the state just by activiting a1,b1, c1 and d1, then the graph is appearing withhin few min..But I wnated to avoid those large expression of a1,b1, c1 and d1. Pl suggest me what to do.
Thanks.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改，无法完成操作。请重新加载页面以查看其更新后的状态。
Translated by 
您也可以从以下列表中选择网站：
如何获得最佳网站性能
选择中国网站（中文或英文）以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)