Is My Code Correct for Computing This Double Sum?

[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".

    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

 

[DrClaude]

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.

$$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##).

 

[member 428835]

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))##.

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))##?

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.

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?

 

[DrClaude]

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.

##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.

##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.

Thanks for replying! So what do you think now?

As far as I can tell, the code you wrote corresponds to your equation.

 

[member 428835]

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!

 

[DrClaude]

Dang it! I'm trying to reproduce results from a paper

Reference?

 

[member 428835]

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!

 

[DrClaude]

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

 

[member 428835]

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

Thanks a lot for your time and advice!