script - matlab integer

Obtener valores integrados de vectores

Question

Estoy tratando de integrar la función seno. Mi objetivo es obtener no solo el valor del área entre una cierta distancia sino los valores específicos del curso integrado.

Una forma de lograr esto es mediante el uso de cumtrapz. Quiero obtener el mismo resultado usando integral o cuádruple. Entonces, me pregunto si hay algo así como "cumquad".

Traté de escribir algo para mí, pero funciona muy lento y parece ser incluso peor que cumtrapz. Más tarde quiero integrar datos medidos. Entonces no será tan simple como un seno.

Aquí está mi código actual:

a = 0; b = 10;
x = a:0.1:b;
y = 2*sin(3*x);
pp = spline(x,y);
y2=zeros(1,length(y));
y3=zeros(1,length(y));
y2(1)=integral(@(x)ppval(pp,x),x(1),x(2));
y3(1)=integral(@(x)ppval(pp,x),x(1),x(2));
for a=2:(length(y)-1)
y2(a) = y2(a-1)+integral(@(x)ppval(pp,x),x(a-1),x(a));
y3(a) = y3(a-1)+quad(@(x)ppval(pp,x),x(a-1),x(a));
end
y4=cumtrapz(x,y);
% y5=cumsum(y);
plot(x,y)
hold on
plot(x,y2,''-ro'')
plot(x,y3,''-kx'')
plot(x,y4,''g'')
syms x % compare with analytical result
ya=2*sin(3*x);
ya5=int(ya)+(2/3);
ezplot(x,ya5)

Answer 1

Usando integral

No creo que haya una forma de que MATLAB devuelva la integral a lo largo de la ruta, por lo que está en lo correcto al realizar la integración un Δx a la vez. La lentitud proviene del ciclo y el siguiente reinicio de cada llamada integral . Puede evitar el ciclo presentando la integral en cada intervalo como una función de valor vectorial.

Las matemáticas

Supongamos que dividimos x en intervalos N-1 con N límites totales y denotamos un límite de intervalo como x _n donde n & in; {1,2,3 ..., N} tal que x ₁ ≤ x ₂ ≤ x ₃ ... ≤ x _N. Entonces cualquier integral sobre el intervalo sería

Usando la u -sustitución:

La integral se convierte en:

donde Δx _n = x _n - x _n-1

El código

Entonces, ahora, podemos plantear la integración de intervalo de cualquier función especificando el límite inferior x _n-1 , especificando el ancho de intervalo Δx , e integrando de 0 a 1 . La mejor parte es que si el límite inferior y el ancho de intervalo son vectores, podemos crear una función con valor vectorial en términos de uy tener una integración integral con la opción ''ArrayValued'' = true .

x = a:0.1:b;
xnm1 = x(1:end-1);
dx = x(2:end) - xnm1;
fx = @(x) 2*sin(3*x);
f = @(u) dx .* fx(dx*u+xnm1);
y = cumsum([0,integral(@(u)f(u),0,1,''ArrayValued'',true)]);

El cumsum explica el hecho de que cada integral en un intervalo dado necesita tener el valor del intervalo anterior agregado a él.

En mi máquina, esto es por lo menos un orden de magnitud más rápido que la versión en bucle y mejora a medida que aumenta el conteo de intervalos.

Usando ode45

El uso también puede usar ode45 para realizar la integración. No es tan eficiente como el método integral , pero puede ser más fácil conceptualmente y verse más limpio. De hecho, ode45 es aproximadamente 10 veces más lento que el método integral anterior cuando se requiere para devolver un error absoluto a la par con el de integral .

a = 0;
b = 10;
% These options are necessary to approach the accuracy of integral
opt = odeset(''RelTol'',100*eps(),''AbsTol'',eps());
sol = ode45(@(x,y) 2*sin(3*x),[a,b],0,opt);
x = a:0.01:b;
yint = deval(sol,x);