How do you solve for a double sum

조회 수: 5 (최근 30일)
nine Yaron
nine Yaron 2019년 11월 13일
답변: jana nassereddine 2021년 4월 29일
Dear everyone:
I was obsessed with a problem for a long time,here is my equation and date.txt.
My calculation result is quite different from the result in the literature. There must be a problem in my code but I don't know what happened. I'm a newbie, somebody please help me. thanks a lot advance!
clear;
format long g
z=load('d:\20191106\a\a1.txt');
n=length(z);
h=0;
i=0.0195314;
m=100
for m=1:512
for j=1:n
for i=1:n-m
h=h+(z(i+m,j)-z(i,j))^2;
end
end
h=1/(n*(n-m))*h;
rr(m)=log(m*i);
hh(m)=log(h);
end
plot(rr,hh,'bo')
My calculations
literatures
  댓글 수: 2
Shubham Gupta
Shubham Gupta 2019년 11월 13일
You have only shared one equation but you are doing lot of other process for e.g. at the end of the 3rd loop about which you have not shared any information. So, I am not sure where should I look for the problem?
It maybe a problem with 'z', it maybe problem with defining 'm' & redefining it inside the loop, it maybe a problem the process at the end of 3rd loop. Possibilies are lot but there is no way I can check that since there is no information about it.
nine Yaron
nine Yaron 2019년 11월 13일
Thank you very much for your comments! I am sorry for lack of information. This equation H(r) is the height-height correlation function of a 3d AFM image, where z(i,j) is the height of the pixel(i,j); m is the number of pixel used for calculation; l is the horizontal length between two pixel. i need the relationship of H(r) and r=ml.

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

채택된 답변

Fabio Freschi
Fabio Freschi 2019년 11월 13일
There are some things that are vgue, like the x-axis, in any case, this code shows something similar to the desired result
clear variables, close all;
z=load('a1.txt');
N = size(z,2);
mMax = size(z,1)-1;
H = zeros(mMax,1);
for m = 1:mMax
for j = 1:N
for i = 1:N-m
H(m) = H(m)+(z(i+m,j)-z(i,j))^2;
end
end
H(m) = 1/(N*(N-m))*H(m);
end
figure,semilogx(1:mMax,log10(H),'bo')
  댓글 수: 1
nine Yaron
nine Yaron 2019년 11월 13일
Oh My Goodness!!! amazing!! it worked! Thank you so much!!! I am really appreciate it for your help! Have a good day!

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

추가 답변 (4개)

Steven Lord
Steven Lord 2019년 11월 13일
h=0;
i=0.0195314;
m=100
for m=1:512
for j=1:n
for i=1:n-m
h=h+(z(i+m,j)-z(i,j))^2;
end
end
h=1/(n*(n-m))*h;
rr(m)=log(m*i);
hh(m)=log(h);
end
Let's take a look at this section of your code. When m becomes 2, what value should h contain? What value does h contain?
You want to put the "h = 0;" line inside your loop over m.
You also want to handle the case where n is equal to m. In that case, your loop over i executes no times and i is []. This causes problems on your line where you assign to rr(m). You probably want to use (n-m) [the last value i would take when the loop executed, which it doesn't when n is equal to m.]
But anyway, we can probably simplify that code a bit. Since the limits of the loop over i don't depend on j at all, we can swap the order of the loops.
for i=1:n-m
for j=1:n
h=h+(z(i+m,j)-z(i,j))^2;
end
end
If we use array operations and sum rather than operating on each scalar element of z independently and using plus, we can eliminate the innermost loop entirely. We do need to be a little careful, since you used length on a matrix and that may not do what you expect if z is taller than it is wide. Instead you should use size to get the number of rows and columns.
for i=1:n-m
v = z(i+m, :)-z(i, :);
h = h + sum(v.^2);
end
You probably can even eliminate that loop if you think of using array operations on submatrices rather than on vectors, but I'll leave that as an exercise for the reader.
  댓글 수: 1
Fabio Freschi
Fabio Freschi 2019년 11월 13일
Yes, I wanted to keep the structure of the code closer to that of the OP. I thought he had more chances to understand. But your code is definitely better from a Matlab style point of view

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


jana nassereddine
jana nassereddine 2021년 4월 29일
hello, can anyone help me solve this,

jana nassereddine
jana nassereddine 2021년 4월 29일
can someone tell me why this code is not working?
clear, clc
f=0; G=0 ;
x= [0; 4]; y= [0; 2]; H=[0; 2]; E=[0; 3];
z=h;
%perpendiculaire surfaces
for l=1:2
for k=1:2
for j=1:2
for i=1:2
Term1= (y(j)-H(k))*(x(i)^2+E(l)^2)^0.5 ;
Term2= atan((y(j)-H(k))/(x(i)^2+E(l)^2)^0.5);
Term3= -0.25*(x(i)^2+E(l)^2-(y(j)-H(k))^2) ;
Term4= ln(x(i)^2+E(l)^2+(y(j)-H(k))^2);
G= (1/(2*pi))*(Term1*Term2+Term3*Term4);
f=f + (-1)^(i+j+k+l)*G;
end
end
end
end
f12 = (1/(x(2)-x(1))*(y(2)-y(1)))*f;

jana nassereddine
jana nassereddine 2021년 4월 29일
pleaseeeee

카테고리

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

태그

제품


릴리스

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by