Numerical Integration Involving Elliptical Functions
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Nicholas Davis
am 16 Feb. 2021
Kommentiert: Alan Stevens
am 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 Kommentare
Akzeptierte Antwort
Alan Stevens
am 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 Kommentare
Alan Stevens
am 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
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Calculus finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!