A 4D Nested numerical integration

2 次查看（过去 30 天）
Taqu3
Taqu3 2022-4-22
评论： Torsten 2022-4-22
Hello. I am new to matlab. I am trying to integrate a function over time and 3-space coordinate using vpaintegral command given below by the S12int1 command. I have waited tens of hours but did not get any result. When it is over 2D given by S12int2 I get the answer in few seconds. In 3D I get a result in about 10 minutes. I list my questions below. Any help is very much appreciated.
  1. What is the stragegy to evaluate multidimensional integrals in general ? Do you recommend vpaintegral over integrate command
  2. In vpaintegral command I can set the accuracy by setting for instance 'RelTol',1e-15,'AbsTol',1e-20. What I observed is if the integrand's magnitude is much larger than the value set by the tolerance values above, integration takes forever. Why is that so ? Do I need to adjust tolerance values in accordance with the integrand's magnitude ?
  3. When using nested vpaintegral commands Do I need to specify error tolerances for each vpaintegral command ? Or just for the inner most integral ? A related question is do I need to put the command 'ArrayValued',true in each nest ?
  4. Is there any other approach ? I am condering turning integrals into summations and perform them in nested for loops. Do you think this would increase the error ?
Thank you!
close all;
%clc;
clear all;
syms t
syms phi
syms z
syms r
tic
zr=1/100; omega=1/100; tau=1; c=1; E0=1/10; kg=2; phik=pi/4; thetak=pi/7;
S12(t,r,z,phi) = exp(-1i*c*kg*t+1i*kg*z*cos(thetak)+ 1i*kg*r*cos(phi-phik)*sin(thetak)) ...
*exp(-(t+z/c)^2/tau^2)*exp(-2*(t-z/c)^2/tau^2)*zr^3/(z^2+zr^2)^(3/2)*exp(-3*r^2*zr*omega/(2*c*(z^2+zr^2)))...
*cos((z/c-t)*omega+r^2*z*omega/(2*c*(z^2+zr^2))-atan(z/zr))*cos((z/c-t)*omega+r^2*z*omega/(2*c*(z^2+zr^2))-atan(z/zr)) ...
*cos((-z/c-t)*omega-r^2*z*omega/(2*c*(z^2+zr^2))+atan(z/zr))*r;
%S12int1= vpaintegral(vpaintegral(vpaintegral(vpaintegral(S12(t,r,z,phi),t,[-3,3]),z,[-10,10]),r,[0,5]),phi,[0,2*pi]);
S12int2= vpaintegral(vpaintegral(vpaintegral(S12(t,r,z,pi/3),t,[-3,3]),r,[0,10]),z,[-10,10],'ArrayValued',true);
S12int2
%double(S12(1,1,1,1))
timeElapsed = toc
  5 个评论
显示 3更早的评论
Taqu3
Taqu3 2022-4-22
Is there a tutorial on how to use the file you have linked above ? Btw, I have tried the following
close all;
%clc;
clear all;
zr=1/100; omega=1/100; tau=1; c=1; E0=1/10; kg=2; phik=pi/4; thetak=pi/7;
S12 = @(t,r,z,phi) exp(-1i.*c.*kg.*t+1i.*kg.*z.*cos(thetak)+ 1i.*kg.*r.*cos(phi-phik).*sin(thetak)) ...
.*exp(-(t+z/c).^2/tau.^2).*exp(-2.*(t-z/c).^2/tau.^2).*zr.^3/(z.^2+zr.^2).^(3/2).*exp(-3.*r.^2.*zr.*omega/(2.*c.*(z.^2+zr.^2)))...
.*cos((z/c-t).*omega+r.^2.*z.*omega/(2.*c.*(z.^2+zr.^2))-atan(z/zr)).*cos((z/c-t).*omega+r.^2.*z.*omega/(2.*c.*(z^.2+zr.^2))-atan(z/zr)) ...
.*cos((-z/c-t).*omega-r.^2.*z.*omega/(2.*c.*(z.^2+zr.^2))+atan(z/zr)).*r;
%integral(@(phi)arrayfun(@(phi)integral3(@(t,r,z)S12(t,r,z,phi),-10,10,0,5,-10,10),phi),0.01,3,'AbsTol',1e-20);
integral(@(phi)integral3(@(t,r,z)S12(t,r,z,phi),-10,10,0,5,-10,10),0,3,'ArrayValued',true,'AbsTol',1e-20)
Both of them gave the warning messages such as
Warning: Matrix is singular to working precision.
Arrays have incompatible sizes for this operation.
Torsten
Torsten 2022-4-22
If you open integralN, you'll see some comments and examples from line 1 to line 75.
Start with simpler examples than yours to see how the code works.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Function Creation 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by