Numerical triple integration of multivariate equations

225 Views Asked by At

The form of my problem is as follows:

$\Psi=C\int_{-\infty}^\infty\int_{-\infty}^\infty\left(f(x,y)\times \left(\int_0^t g(x,t)dt\right)dxdy\right)$

I have already attempted numerical solutions using Matlab's integrate3, without success.

In any case I doubt integrate3 is the right approach as that is meant for problems of the form:

$\Psi=\int_a^b\int_c^d\int_e^f f(x,y,z)dxdydz$

Similarly, attempts made with scipy's integration toolkit have also not borne fruit.

I have attempted to first calculate $\int_0^tg(x,t)$ at discrete $x$ values and then place it in $\Psi$ but for rather obvious reasons that does not work either.

Additionally, $g(x,t)$ cannot be factored in the form of $h(x)\times i(t)$, which might have allowed for a by parts solution which might be integrated symbolically (for x) and numerically for t.

Also $\int_0^t\int_{-\infty}^{\infty} g(x,t)$ has singularities at multiple points.

Is there a canonical method of solving this?

UPDATE

The equations being modeled are given in the paper by Shen et. al. (in spherical coordinates, hence the reduction of integrals) as:

$\Theta=\frac{\theta}{t_c}\int_0^t\frac{1}{1+\frac{2t^\prime}{t_c}}\left[1-\mathrm{exp}\left(\frac{-2mg}{1+2\frac{t^\prime}{t_c}}\right)dt^\prime\right]$

and

$U_p=C\int_0^\infty \mathrm{exp}[-(1+\mathrm{i}V)g]\mathrm{exp}(-\mathrm{i}\Theta) dg$

With the approximation $\mathrm{exp}(-\mathrm{i}\Theta)\approx1-\mathrm{i}\Theta$

$U=C\int_0^{\infty}(1-\mathrm{i}\Theta)\mathrm{exp}[-(1+\mathrm{i}V)g]dg$

The form I've explained in the original question is obtained on converting to Cartesian coordinates.

UPDATE 3

Using symbolic math I get

for t=1:1000:Tmax
phi = (theta/tc)* int((1/(1+2*(tp/tc)))*[1-exp((2*m*g)/(1+(2*            
(tp/tc))))],tp,0,t);
phi=subs(phi)
preU=exp(-(1+1i*V)*g)*exp(-1i*phi);

Urun=int(preU,g,0,inf); % Evaluated at t=0

tt(t+1)=t+1;
end

Io= 3/193 - 84i/693;

Ir=abs(Io);

Uplot = (vpa(abs(Urun))./Ir).^2

However this returns a value of Uplot in the form abs(numeric::int()) which cannot be plotted.

Additionally, attempting to numerically integrate the second portion of the integral (the portion returning the abs(numeric::int()) form) I get an array of NaN's.

Needless to say, still no plot.

for t=1:1:Tmax
phi = (theta/tc)* int((1/(1+2*(tp/tc)))*[1-exp((2*m*g)/(1+(2*
(tp/tc))))],tp,0,t);
phi=subs(phi);
preU=exp(-(1+1i*V)*g)*(1-1i*phi);

UrunF=matlabFunction(preU)
UrunNum(t)=integral(UrunF,0,inf);

tt(t+1)=t+1;
end

Any suggestions?

1

There are 1 best solutions below

0
On

If I understand your problem correctly, to compute $\int_{-\infty}^\infty\int_{-\infty}^\infty\left(f(x,y)\times \left(\int_0^t g(x,\tau)d\tau\right)dxdy\right), $ I would start by defining a subfunction that computes $G(x, t) \equiv \int_0^t g(x,\tau)d\tau, $ and then, for each $t,$ integrate $f(x,y)\times G(x, t)$ over the plane with integral2.

As an example with some arbitrary functions $f$ and $g:$

f = @(x, y) exp(-(x.^2 + y.^2));
g = @(t, x) abs(x) .* exp(-t);
t_values = 0:10;
n = numel (t_values);
S = nan (n, 1);
for i = 1:n
  t = t_values(i);  
  G = @(x) arrayfun(@(x) integral(@(tau) g(tau, x), 0, t), x);
  S(i) = integral2 (@(x, y) f(x, y) .* G(x), -Inf, Inf, -Inf, Inf);
end
plot (t_values, S)