how plot amplitude of mode from eigenvector?

4 Ansichten (letzte 30 Tage)

Habib am 22 Nov. 2015

0
Verknüpfen

Direkter Link zu dieser Frage

https://de.mathworks.com/matlabcentral/answers/256896-how-plot-amplitude-of-mode-from-eigenvector

Hi every body.

For solving an integral (this integral is Huygens Fresnel integral and it is useful for find modes of a laser resonator), I used Gauss-Legendre quadrature method and in final I got its eigenvectors and eigenvalues. but I couldn't plot the mode's amplitude and phase. If anyone has a experience about it, I hope help?

these codes are what is that I wrote.

 clc
clear all
close all
a=0;
b=1;
N=200;
N=N-1;
N1=N+1; N2=N+2;
xu=linspace(a,b,N1)';
% Initial guess
y=cos((2*(0:N)'+1)*pi/(2*N+2))+(0.27/N1)*sin(pi*xu*N/N2);
% Legendre-Gauss Vandermonde Matrix
L=zeros(N1,N2);
% Derivative of LGVM
Lp=zeros(N1,N2);
% Compute the zeros of the N+1 Legendre Polynomial
% using the recursion relation and the Newton-Raphson method
y0=2;
% Iterate until new points are uniformly within epsilon of old points
while max(abs(y-y0))>eps
    L(:,1)=1;
    Lp(:,1)=0;
    L(:,2)=y;
    Lp(:,2)=1;
    for k=2:N1
        L(:,k+1)=( (2*k-1)*y.*L(:,k)-(k-1)*L(:,k-1) )/k;
    end
    Lp=(N2)*( L(:,N1)-y.*L(:,N2) )./(1-y.^2);
    y0=y;
    y=y0-L(:,N2)./Lp;
end
% Linear map from[-1,1] to [a,b]
x=(a*(1-y)+b*(1+y))/2;
% Compute the weights
w=(b-a)./((1-y.^2).*Lp.^2)*(N2/N1)^2;
LL=40.92;
A=1; B=LL; C=0;
n=N+1; m=n; s=n;
landa=10.6e-4;
k=2*pi/landa;
r1=linspace(a,b,n);
r2=linspace(a,b,m);
kernel1=zeros(n,m);
kernel2=zeros(n,m);
H1=zeros(n,m);
H2=zeros(n,m);
T1=zeros(1,m);
T2=zeros(1,m);
R=100*LL;
teta0=12.22e-3*pi/360;
T1(:)=exp(-2*1i*k*teta0*r1(:));
T2(:)=exp(-1i*k*r1(:).^2./R);
for s=1:n
    kernel1(:,s)=-1i*(k/B)*r2(s).*besselj(0,r2(s).*r1(:)*k/B).*exp(1i*k*(A*r2(s).^2+B*r1(:).^2)/B/2);
   H1(:,s)=kernel1(:,s).*T1(:);
end
for s=1:m
    kernel2(:,s)=-1i*(k/B)*r1(s).*besselj(0,r1(s).*r2(:)*k/B).*exp(1i*k*(A*r1(s).^2+B*r2(:).^2)/B/2);
    H2(:,s)=kernel2(:,s).*T2(:);
end
W=diag(w);
P1=H1*W;
P2=H2*W;
HH=P1*P2;
[V,D]=eig(HH);