Matrix inversion and summation

조회 수: 3 (최근 30일)
BeeTiaw
BeeTiaw 2018년 6월 29일
댓글: BeeTiaw 2018년 6월 30일
Hi all, can anyone help me to explain the mistake in my code below? The problem is as follow:
  • The following parameter is assumed,
lma = 3;
ma = zeros(1,lma);
ma(2) = -0.2;
ma(lma) = -0.01;
  • The parameter is then used to calculate the coordinate z using the following equation,
where z = x + i*y (complex number)
dumZ = 0;
sumdumZ = 0;
theta = linspace(0,2*pi,37);
xi = exp(1i.*theta);
for ii=1:length(theta)
for jj=1:lma
dumZ = ma(jj).*xi(ii).^jj;
sumdumZ = sumdumZ + dumZ;
end
zcon(ii) = 1./xi(ii)+sumdumZ;
dumZ = 0;
sumdumZ = 0;
end
zcon=zcon';
The following profile is obtained by plotting zcon,
  • Now, using the hypothetical data above, i.e. the zcon, I want to back-calculate the variable "ma" or "alpha" as shown in the equation above. To do this, I first created a matrix A, which essentially a matrix containing all the constants. The following code is used to generate matrix A with 370 terms,
lx = length(zcon);
N = 10*37;
A = zeros(lx,N);
for ii=1:lx
for jj=1:N
A(ii,jj) = xi(ii).^jj;
end
end
A=[1./xi' A]; % insert a new colum to the left containing the 1/xi(ii) into the existing matrix A;
  • Then the variable "ma" or "alpha" can be easily obtained using the following operation in Matlab,
maiter = A\zcon;
which give me a [370 x 1] matrix. To validate this, I again calculate the "predicted" z using the following code:
zdirect = A*maiter;
Which results an exactly the same curve (as expected) as shown below,
Now, coming to my question. Instead of using a direct calculation as shown above, i.e. the code "zdirect = A*maiter, instead I would like to manually use summation to calculate the "z" for each "theta" using the code below. Somehow, I could not get it right and cannot find the mistake. Below is the code that I used to calculate "z"
dumZ = 0;
sumdumZ = 0;
pow = [-1 1:1:N];
for ii=1:length(theta)
for jj=1:length(maiter)
dumZ = maiter(jj).*xi(ii).^pow(jj);
sumdumZ = sumdumZ + dumZ;
end
ziter(ii) = sumdumZ;
dumZ = 0;
sumdumZ = 0;
end
Instead, this is what I got,
Can someone please help me?

채택된 답변

David Goodmanson
David Goodmanson 2018년 6월 29일
Hi Bee,
Are you plotting the right variable? I ran your code, and zcon, zdirect and ziter all reproduce the same curve. Only the simple exponential variable xi gives the red circle in the second plot.
  댓글 수: 4
David Goodmanson
David Goodmanson 2018년 6월 30일
Hi Bee,
Below are the first few rows of [zcon zdirect ziter]. The only change I made to your code was
zcon = zcon.' instead of zcon = zcon'
so that the complex conjugate of zcon is not taken. That means that the table I get is the complex conjugate of yours. I also added
ziter = ziter.'
so that ziter is a column vector. Of course neither of these changes/additions affects the shape of the plot, 'triangle' vs. circle.
0.7900 + 0.0000i 0.7900 + 0.0000i 0.7900 + 0.0000i
0.7882 - 0.2471i 0.7882 - 0.2471i 0.7882 - 0.2471i
0.7815 - 0.4792i 0.7815 - 0.4792i 0.7815 - 0.4792i
0.7660 - 0.6832i 0.7660 - 0.6832i 0.7660 - 0.6832i
0.7363 - 0.8484i 0.7363 - 0.8484i 0.7363 - 0.8484i
0.6862 - 0.9680i 0.6862 - 0.9680i 0.6862 - 0.9680i
0.6100 - 1.0392i 0.6100 - 1.0392i 0.6100 - 1.0392i
0.5039 - 1.0633i 0.5039 - 1.0633i 0.5039 - 1.0633i
0.3666 - 1.0446i 0.3666 - 1.0446i 0.3666 - 1.0446i
0.2000 - 0.9900i 0.2000 - 0.9900i 0.2000 - 0.9900i
0.0093 - 0.9077i 0.0093 - 0.9077i 0.0093 - 0.9077i
-0.1975 - 0.8061i -0.1975 - 0.8061i -0.1975 - 0.8061i
-0.4100 - 0.6928i -0.4100 - 0.6928i -0.4100 - 0.6928i
-0.6167 - 0.5741i -0.6167 - 0.5741i -0.6167 - 0.5741i
-0.8058 - 0.4545i -0.8058 - 0.4545i -0.8058 - 0.4545i
BeeTiaw
BeeTiaw 2018년 6월 30일
Hi David Goodmanson,
Thanks for your answer!
I got it correct now after change it as per your suggestion. Thanks a lot!
It seems that it has something to do with consistency of the conjugate/non-conjugate of the complex number. Although I still do not clearly understand it.
But I would consider the problem is solved now.

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Mathematics에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by