Trapezoidal Rule Involving Elliptical Integrals/Functions

Hi. I'm trying to write a code to integrate a function via the trapezoidal rule. I can't seem to get my graph to produce the correct output. I have attached an image of the expected output. The function I am attempting to integrate is dphi = alpha/p. I am also unsure if I used the correct formula for the trapezoidal rule, so any clarity there would be appreciated. The first figure is exactly what I am looking for but the second figure cannot produce the correct result. I don't assume any fault in the math here as plotting phi versus x should be simple. I have attempted this issue before with Matlab's integral command but it still did not produce the correct plot.

clear all;

% Limits of Integration
a = 0;   b = 1;   n = 100;
h = (b-a)/n;

% Prerequisites
x = linspace(0,1,n);
m = 0.999999129426574; 
J = 1;
L = 0.04;
[K,E] = ellipke(m);
r = zeros(1,numel(x));
dphi = zeros(1,numel(x));
trapphi = zeros(1,numel(x));

% Function
s = 1 + 8*J^2*L^2*E*K;
t = 8*J^2*L^2*K^2;
alpha = (1/sqrt(2)*L)*sqrt(s*(s-(1-m)*t)*(s-t));

for i = 1:n-1
    u = 2*J*K*x(i);
    [SN] = ellipj(u,m);
    r(i) = s - t + t * m * SN^2;
    dphi(i) = alpha/r(i);   % Function to integrate
    trapphi(i) = (sum(dphi) - (dphi(i)+dphi(i+1))/2)*h;
    phi(i) = mod(trapphi,2*pi);
end

% Plotting
figure(1)   % phi versus x
plot(x,phi,'-k')
xlim([0,1]);
ylim([0,1]);
figure(2)
plot(x,r,'-k')   % r versus x
xlim([0,1]);

[Allen]

Nicholas,

Try looking at documentation on MATLAB's built-in trapezoidal functions (trapz and cumtrapz).

[Mathieu NOE]

I could get the shape but not the range for phi. I still need to find out

clear all;
% Limits of Integration
a = 0;   b = 1;   n = 100;
h = (b-a)/n;
% Prerequisites
x = linspace(0,1,n);
m = 0.999999129426574; 
J = 1;
L = 0.04;
[K,E] = ellipke(m);
r = zeros(1,numel(x));
dphi = zeros(1,numel(x));
trapphi = zeros(1,numel(x));
% Function
s = 1 + 8*J^2*L^2*E*K;
t = 8*J^2*L^2*K^2;
alpha = (1/sqrt(2)*L)*sqrt(s*(s-(1-m)*t)*(s-t));
for i = 1:n % here
    u = 2*J*K*x(i);
    [SN] = ellipj(u,m);
    r(i) = s - t + t * m * SN^2;
    dphi(i) = alpha/r(i);   % Function to integrate
%     trapphi(i) = (sum(dphi) - (dphi(i)+dphi(i+1))/2)*h;
%     phi(i) = mod(trapphi(i),2*pi); % here
end
phi = cumtrapz(x,dphi);

% Plotting
figure(1)   % phi versus x
plot(x,phi/(2*pi),'-k')
% xlim([0,1]);
% ylim([0,1]);
figure(2)
plot(x,r,'-k')   % r versus x
xlim([0,1]);

[David Goodmanson]

Hi Nicholas,

the top figure looks good as you say, but from the shape of the bottom figure, it does not appear to be the integral of (alpha / r), even after correcting for the size disagreement by multiplying by an appropriate constant . Could you provide a source and some more background for these equations?

[Nicholas Davis]

Hi David,

The code Mathieu wrote is more than sufficient for my needs, but I was hoping I could fit the range perfectly! I have attached the sources for my equations. The equation of (alpha/r) is on page 6, and it's labeled (41). The figures I'm trying to replicate are from figure 4, a and c, on page 7. You will find the highlights very helpful. My notation is slightly different from those on the attached paper, so be aware of that. Thanks!

[David Goodmanson]

Hi Nicholas,

I took a look at the paper. Things agree pretty well except that in the expression for alpha you have

(1/sqrt(2)*L)     rather than (1/(sqrt(2)*L))

Making that change brings the answer up by a factor of 625, which is at least into the realm of possibility. The result for phi is about 10.3 which is a bit too large, and the shape still does not look like fig 4b.

You can streamline the code by using

u = 2*J*K*x;
sn = ellipj(u,m);
r = s - t + t * m * sn.^2;
dphi = alpha./r;   % Function to integrate

phi = cumtrapz(x,dphi);

in which case you can drop the preallocation of several variables.

[Mathieu NOE]

hello guys

reading the publication was no easy task for me - a bit far from my usual area of expertise....

could not really figyre out where there could be a bug

nevertheless, I just noticed that the figures in the publication seem to differ a bit in terms of initial and final slope

the publication plot seem to show more vertical and straight lines vs what this code generates ( with above corrections) and also I increased the number of points from 100 to 1000 to see if refining the grid would improve the results (not that much indeed

so my believe is that the "r" function has not the right slope somewhere - but I ignore the reasons why (my 2 cents)

clear all;
% Limits of Integration
a = 0;   b = 1;   n = 1000;
h = (b-a)/n;
% Prerequisites
x = linspace(a,b,n);
m = 0.999999129426574; 
J = 1;
L = 0.04;
[K,E] = ellipke(m);
r = zeros(1,numel(x));
dphi = zeros(1,numel(x));
trapphi = zeros(1,numel(x));
% Function
s = 1 + 8*J^2*L^2*E*K;
t = 8*J^2*L^2*K^2;
alpha = (1/(sqrt(2)*L))*sqrt(s*(s-(1-m)*t)*(s-t));

u = 2*J*K*x;
sn = ellipj(u,m);
r = s - t + t * m * sn.^2;
dphi = alpha./r;   % Function to integrate

phi = cumtrapz(x,dphi);
% Plotting
figure(1)   % phi versus x
plot(x,phi/(2*pi),'-k')
% xlim([0,1]);
% ylim([0,1]);
figure(2)
plot(x,r,'-k')   % r versus x
xlim([0,1]);

[Bjorn Gustavsson]

Maybe this is related to your discrepancies: I recall having problems using cumtrapz on functions where the limit goes to infinity but the integral converges - this was also for denominators of the type sqrt(a^2-x.^2).

[Nicholas Davis]

Hi all, I was able to get the plots from the paper perfectly with everyone's help. My professor originally advised I use the mod function to get my phi, but that was actually bad advice. In the new code I simply adjusted phi so instead of using the mod function, it was phi/(2*pi) and that gave me the correct range. For those interested, I was able to get all of the plots in figure 4 from the paper more or less perfectly. I need to develop a better newton's method code for finding the proper m values and thus getting better results, but I digress. I have attached the final code for all of you to marvel at what you helped create. Thank you everyone.

clear all
format long
% Constants ---------------------------------------------------------------
n = 1000;
x = linspace(0,1,n);
J = input('Enter J value: ');
L = 0.04;
r = zeros(1,numel(x));
dphi = zeros(1,numel(x));

% Trapz Command Integration -----------------------------------------------
if J == 1
    for i = 1:n
        m = 0.999999129426574;
        [K,E] = ellipke(m);
        u = 2*J*K*x(i);
        [SN] = ellipj(u,m);
        s = 1 + 8*J^2*L^2*E*K;
        t = 8*J^2*L^2*K^2;
        r(i) = s - t + t * m * SN^2;
        alpha = (1/(sqrt(2)*L))*sqrt(s*(s-(1-m)*t)*(s-t));
        dphi(i) = alpha/r(i);
    end
trapzphi = cumtrapz(x,dphi);
phi = mod(trapzphi,2*pi);

elseif J == 2
    for i = 1:n
        m = 0.997999129426574; % 0.993524799088048;
        [K,E] = ellipke(m);
        u = 2*J*K*x(i);
        [SN] = ellipj(u,m);
        s = 1 + 8*J^2*L^2*E*K;
        t = 8*J^2*L^2*K^2;
        r(i) = s - t + t * m * SN^2;
        alpha = (1/(sqrt(2)*L))*sqrt(s*(s-(1-m)*t)*(s-t));
        ialpha = abs(alpha);
        dphi(i) = ialpha/r(i);
    end
    trapzphi = cumtrapz(x,dphi);
    phi = trapzphi/(2*pi); % mod(trapzphi,2*pi);    
end

% Plotting ----------------------------------------------------------------
clf
figure(1)   % r versus x
plot(x,r,'-k')
figure(2)
plot(x,phi,'-k')   % phi versus x

[Mathieu NOE]

Congrats to the team work !