Numerical Integration Involving Elliptical Functions

Hi.

I'm trying to plot a series of equations. The first of which will be fit into figure 1 (x versus p) and the other being fit into figure 2 (x versus phi). To derive phi based on the paper I was given requires integrating, but I am running to an error that I cannot seem to resolve. I understand the error, but I can't seem to find the change that would fix it. I have attached the code. I am also open to any other changes to better the code. Thanks. This is the error:

Error using integralCalc/finalInputChecks (line 526)

Output of the function must be the same size as the input. If FUN is an array-valued integrand, set the 'ArrayValued' option to true.

Error in integralCalc/iterateScalarValued (line 315)

finalInputChecks(x,fx);

Error in integralCalc/vadapt (line 132)

[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)

[q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 88)

Q = integralCalc(fun,a,b,opstruct);

Error in true_carr_fig4 (line 20)

phi = integral(fun,0,x);

clear all;
% Initial Values
m = 0.999999119426574;   % Value from Newton's Method
scale = 0.04;   J = 1;
[K,E] = ellipke(m);

% Workspace 4.0
for x = 0:0.01:1
    u = 2*J.*K.*x;
    [SN] = ellipj(u,m);
    JSN = SN.^2;
    s = 1 + 8 * J^2 * scale^2 .* E .*K;
    t = 8 * J^2 * scale^2 .* K.^2;
    p = s - t + t .* m .* JSN;
    alpha = (1/(sqrt(2)*scale)).*sqrt(s.*(s-(1-m).*t).*(s-t));
    figure(1)
    plot(x,p,'.k');
    fun = @(m) alpha./p;
    ph = integral(fun,0,x);
    phi = ph/2*pi;
    figure(2)
    plot(x,phi,'.k');
    hold on
end
hold off

%xlim([-0.05,1.05]);
%ylim([0,1.2]);

Answer 1

[Alan Stevens]

More like this perhaps:

% Initial Values
m = 0.999999119426574;   % Value from Newton's Method
scale = 0.04;   J = 1;
[K,E] = ellipke(m);
x = 0:0.01:1;
p = zeros(1,numel(x));
phi = zeros(1,numel(x));
% Workspace 4.0
for i = 1:numel(x)
    u = 2*J.*K.*x(i);
    [SN] = ellipj(u,m);
    JSN = SN.^2;
    s = 1 + 8 * J^2 * scale^2 .* E .*K;
    t = 8 * J^2 * scale^2 .* K.^2;
    p(i) = s - t + t .* m .* JSN;
    alpha = @(m) (1/(sqrt(2)*scale)).*sqrt(s.*(s-(1-m).*t).*(s-t));
    fun = @(m) alpha(m)./p(i);
    ph = integral(fun,0,x(i));
    phi(i) = ph/2*pi;
end

figure(1)
plot(x,p,'.k'),grid
xlabel('x'),ylabel('p')
figure(2)
plot(x,phi,'.k'), grid
xlabel('x'),ylabel('\phi')

[Nicholas Davis]

Thank you Alan. I didn't think of preallocating arrays!

[Nicholas Davis]

While we're at it, I want to compare Matlab's integral and rsums commands, but it gives me a "too many output arguments" error. I figured that preallocating the arrays would offer some benefit, but it did not work. These are the only additions I made:

% Initial Values
format long
m = 0.999999129426574;   % Value from Newton's Method
scale = 0.04;   J = 1;
[K,E] = ellipke(m);
x = 0:0.01:1;
p = zeros(1,numel(x));
phi = zeros(1,numel(x));
riemann = zeros(1,numel(x)); % Added for rsums
% Workspace 4.0
for i = 1:numel(x)
    u = 2*J.*K.*x(i);
    [SN] = ellipj(u,m);
    JSN = SN.^2;
    s = 1 + 8 * J^2 * scale^2 .* E .*K;
    t = 8 * J^2 * scale^2 .* K.^2;
    p(i) = s - t + t .* m .* JSN;
    alpha = @(m) (1/(sqrt(2)*scale)).*sqrt(s.*(s-(1-m).*t).*(s-t));
    fun = @(m) alpha(m)./p(i);
    ph = integral(fun,0,x(i));
    phi(i) = mod(ph,2*pi)/2*pi; %ph/2*pi;
    riemann(i) = rsums(fun,0,x(i)); % Added for rsums
end

figure(1)
plot(x,p,'.k'),grid
xlabel('x'),ylabel('p')
xlim([0,1]);
figure(2)
plot(x,phi,'.k'), grid
xlabel('x'),ylabel('\phi')
ylim([0,1]);   xlim([0,1]);

[Alan Stevens]

rsums works differently. It just produces a histogram, the number of terms used can then be adjusted interactively.

doc rsums

The following works for a few values of x, but you probably won't want to generate 100 histograms all in one go!

% Initial Values
format long
m = 0.999999129426574;   % Value from Newton's Method
scale = 0.04;   J = 1;

x = [0.1 0.5 1];
for i = 1:numel(x)
[K,E] = ellipke(m);
u = 2*J.*K.*x(i);
[SN] = ellipj(u,m);
JSN = SN.^2;
s = 1 + 8 * J^2 * scale^2 .* E .*K;
t = 8 * J^2 * scale^2 .* K.^2;
p(i) = s - t + t .* m .* JSN;
alpha = @(m) (1/(sqrt(2)*scale)).*sqrt(s.*(s-(1-m).*t).*(s-t));
fun =  @(m) alpha(m)./p(i);
figure
rsums(fun,0,x(i));
end