現在この質問をフォロー中です
- フォローしているコンテンツ フィードに更新が表示されます。
- コミュニケーション基本設定に応じて電子メールを受け取ることができます。
Numerical Integration in matlab
3 ビュー (過去 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 で Characters and Strings についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!エラーが発生しました
ページに変更が加えられたため、アクションを完了できません。ページを再度読み込み、更新された状態を確認してください。
Web サイトの選択
Web サイトを選択すると、翻訳されたコンテンツにアクセスし、地域のイベントやサービスを確認できます。現在の位置情報に基づき、次のサイトの選択を推奨します:
また、以下のリストから Web サイトを選択することもできます。
最適なサイトパフォーマンスの取得方法
中国のサイト (中国語または英語) を選択することで、最適なサイトパフォーマンスが得られます。その他の国の 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 (한국어)