이 질문을 팔로우합니다.
- 팔로우하는 게시물 피드에서 업데이트를 확인할 수 있습니다.
- 정보 수신 기본 설정에 따라 이메일을 받을 수 있습니다.
Numerical Integration in matlab
조회 수: 12 (최근 30일)
이전 댓글 표시
AVM
2020년 1월 22일
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
James Tursa
2020년 1월 22일
Please give us the details on what you are trying to integrate, and over what range.
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!오류 발생
페이지가 변경되었기 때문에 동작을 완료할 수 없습니다. 업데이트된 상태를 보려면 페이지를 다시 불러오십시오.
웹사이트 선택
번역된 콘텐츠를 보고 지역별 이벤트와 혜택을 살펴보려면 웹사이트를 선택하십시오. 현재 계신 지역에 따라 다음 웹사이트를 권장합니다:
또한 다음 목록에서 웹사이트를 선택하실 수도 있습니다.
사이트 성능 최적화 방법
최고의 사이트 성능을 위해 중국 사이트(중국어 또는 영어)를 선택하십시오. 현재 계신 지역에서는 다른 국가의 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 (한국어)