Numerical Integration Involving Elliptical Functions
7 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Nicholas Davis
el 16 de Feb. de 2021
Comentada: Alan Stevens
el 17 de Feb. de 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 comentarios
Respuesta aceptada
Alan Stevens
el 16 de Feb. de 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 comentarios
Alan Stevens
el 17 de Feb. de 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
Más respuestas (0)
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!