Trapezoidal Rule Involving Elliptical Integrals/Functions
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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]);
9 comentarios
Respuestas (0)
Ver también
Categorías
Más información sobre Numerical Integration and Differentiation en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!