Is My Code Correct for Computing This Double Sum?

  • MATLAB
  • Thread starter member 428835
  • Start date
  • Tags
    Sum
  • #1
member 428835
Hi PF!

I'm trying to compute $$S_i(\theta) = \sum_{m=odd}^{N}\frac{64}{\pi^3 m \alpha}\frac{k^2}{4k^2-m^2} \sum_{l=0}^N J_l(\alpha m) \xi_l \left( B_l(\theta)\cos(l\theta)+C_l(\theta)\sin(l\theta) \right)$$ where the ##i## subscript appears since $$B_l = \int_0^{2\pi} \psi_i(\theta) \cos(l\theta)\,d\theta \\ C_l = \int_0^{2\pi} \psi_i(\theta) \sin(l\theta)\,d\theta.$$ Once obtaining ##S_i## I'm trying to compute $$[M]:M_{ij} = \int_0^{2\pi}S_i\psi_{j} \, d\theta.$$ Can you confirm that this code computes that? I'm driving myself crazy here. In the code below I call ##S_i## "Mop".
Code:
    for ii = 1:2*N-1
        for mm = 1:2:N
            for l = 0:N
                J = 2*besseli(l,alpha*mm)/(besseli(l-1,alpha*mm)+besseli(l+1,alpha*mm));
                B = trapz(psi(ii,:).*cos(l.*theta))*dtheta;
                C = trapz(psi(ii,:).*sin(l.*theta))*dtheta;
                Mop1 = Mop1+xi(l+1)*J*(B*cos(l*theta)+C*sin(l*theta));
            end% for l

            Mop(ii,:) = Mop(ii,:) + 64/(pi^3*mm*alpha)*(k^2/(4*k^2-mm^2))*Mop1;           
            Mop1 = zeros(1,size(Mop1,2));% clear previous Mop1 data
        end% for mm

        for jj = 1:2*N-1
            M(ii,jj) = trapz(Mop(ii,:).*psi(jj,:))*dtheta;% eq. (2.24b)
        end% for jj
    end% for ii
 
Physics news on Phys.org
  • #2
i am bothered by the fact that Mop1 is set to zero at the end the loop. This doesn't guarantee that Mop1 is zero upon the first iteration of the loop.

In the calculation of ##J_0##, you have the modified Bessel function with a negative ##\nu##. That is not a problem per se, but is it what you actually want?

It would be helpful to know the value of N and the size of the vector theta.

joshmccraney said:
$$S_i(\theta) = \sum_{m=odd}^{N}\frac{64}{\pi^3 m \alpha}\frac{k^2}{4k^2-m^2} \sum_{l=0}^N J_l(\alpha m) \xi_l \left( B_l(\theta)\cos(l\theta)+C_l(\theta)\sin(l\theta) \right)$$
I guess that writing ##B_l(\theta)## and ##C_l(\theta)## is a typo (both coefficients are independent of ##\theta##).
 
  • #3
DrClaude said:
i am bothered by the fact that Mop1 is set to zero at the end the loop. This doesn't guarantee that Mop1 is zero upon the first iteration of the loop.
Sorry, the actual code I use is much longer so I was trying to be compact but clearly left some things out. Before the for ii loop I preallocate ##Mop1## to zero by ##
Mop1 = zeros(1,length(theta))##.

DrClaude said:
In the calculation of ##J_0##, you have the modified Bessel function with a negative ##\nu##. That is not a problem per se, but is it what you actually want?
##J_l(\alpha m) = I_l(\alpha m)/I'_l(\alpha m)##, where ##I## is the modified Bessell function of the first kind. So isn't it true that ##J_0 =
2I_0(\alpha m)/(I_{-1}(\alpha m)+I_1(\alpha m))##?

DrClaude said:
It would be helpful to know the value of N and the size of the vector theta.
##N=10## and theta is a ##1\times 4284## size vector.

DrClaude said:
I guess that writing ##B_l(\theta)## and ##C_l(\theta)## is a typo (both coefficients are independent of ##\theta##).
Yes, definitely a typo on my part, thanks for pointing this out.

Thanks for replying! So what do you think now?
 
  • #4
joshmccraney said:
Sorry, the actual code I use is much longer so I was trying to be compact but clearly left some things out. Before the for ii loop I preallocate ##Mop1## to zero by ##
Mop1 = zeros(1,length(theta))##.
While this leads to a correct result, I don't like the approach for two reasons. First, it makes the loop harder to figure out. Second, it leads to an unnecessary Mop1 = 0 at the end of the loop.

joshmccraney said:
##J_l(\alpha m) = I_l(\alpha m)/I'_l(\alpha m)##, where ##I## is the modified Bessell function of the first kind. So isn't it true that ##J_0 =
2I_0(\alpha m)/(I_{-1}(\alpha m)+I_1(\alpha m))##?
Ok. I just didn't know if what you wrote corresponded to what you intended.

joshmccraney said:
##N=10## and theta is a ##1\times 4284## size vector.
I wanted to check if you were not using too few points per period.

joshmccraney said:
Thanks for replying! So what do you think now?
As far as I can tell, the code you wrote corresponds to your equation.
 
  • Like
Likes member 428835
  • #5
DrClaude said:
As far as I can tell, the code you wrote corresponds to your equation.
Dang it! I'm trying to reproduce results from a paper and I had the correct answer for a scenario that is identical to this in every way except there was no double sum, which is to say there was no sum over odd indices. Since my results were not accurate, I was hoping I messed up the double sum.

I did try a double sum on something like ##\sum_{i=1}^3 i \sum_{j=1}^3 ij## and it worked, but I thought maybe I messed something else up. I suppose the paper could have a mistake, though this seems unlikely. Either way, thanks so much for your input and for taking the time to help me out!
 
  • #6
joshmccraney said:
Dang it! I'm trying to reproduce results from a paper
Reference?
 
  • #7
DrClaude said:
Reference?
Stability of constrained cylindrical interfaces and the torus lift of Plateau–Rayleigh. Written by Josh Bostwick and Paul Steen. The math side of things really starts at equation (2.16). I can reproduce results for the natural boundary condition, section 2.1.1, but using the same code and changing ##\alpha## and the matrix operator ##M## as mentioned above, I am unable to reproduce results for the pinned boundary, section 2.1.2 (##M## and ##\alpha## are the only two differences in these two scenarios).

I think I found a small error in their computation for results for section 2.1.1: see I can reproduce their results for this natural boundary condition but only if I slightly change their ##M## operator in equation (2.19c), specifically by dividing by ##\pi^2## instead of ##\pi##. I re-derived equation (2.19c) and agree that (2.19c) is correct. However, when re-deriving equation (2.21c) I do not get what they have. I can send you details if you're interested ( have them LaTeX'ed up in a separate PDF)? I'd appreciate all the help you're willing to offer.

Again, thanks so much for your time!
 
  • #8
I have had a quick look at the article, but unfortunately, I do not have much time right now to look at this in more details.

I suggest you contact the authors. Some can be very interested in talking about their work :smile:
 
  • #9
DrClaude said:
I have had a quick look at the article, but unfortunately, I do not have much time right now to look at this in more details.

I suggest you contact the authors. Some can be very interested in talking about their work :smile:
Thanks a lot for your time and advice!
 

Similar threads

Replies
3
Views
2K
Replies
7
Views
2K
Replies
6
Views
6K
Replies
12
Views
918
Replies
2
Views
2K
Replies
3
Views
6K
Back
Top