Numerical Integration Involving Elliptical Functions
5 views (last 30 days)
Show older comments
Nicholas Davis
on 16 Feb 2021
Commented: Alan Stevens
on 17 Feb 2021
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]);
0 Comments
Accepted Answer
Alan Stevens
on 16 Feb 2021
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')
3 Comments
Alan Stevens
on 17 Feb 2021
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
More Answers (0)
See Also
Categories
Find more on Calculus in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!