MATLAB not computing integral of an infinite integral

Question

AD 2023-8-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2011692-matlab-not-computing-integral-of-an-infinite-integral

编辑： Dyuman Joshi 2023-8-28

I am trying to compute the integral of the function. However, MATLAB is unabble to compute it. I havetried numerical integration function 'integral' with no results. Can someone please tell how to proceed with this? Thanks in advance!

P_l=50;
v=0.1;
k=15;
Tm=1;
T0=300;
alpha=3.75*10^(-6);
E= 190*10^9;
nu=0.3;
syms x z x_prime z_prime t dT_dx dT_dz;
xm=x-x_prime;
zp=z+z_prime;
zm=z-z_prime;
% Define the terms
T=(P_l*exp(-((sqrt((x-v*t)^2 + (z)^2) + (x-v*t)))/(2*alpha)))/(4*3.14*k*sqrt((x-v*t)^2 + (z)^2)) +T0;
dT_dx=diff(T,x);
dT_dx_prime=subs(dT_dx,[x,z],[x_prime,z_prime]);
dT_dz=diff(T,z);
dT_dz_prime=subs(dT_dz,[x,z],[x_prime,z_prime]);
Gxh=@(xm,zm,zp) (1/(4*pi))*(3*(xm/(xm^2 + zp^2)) + 2*(xm*zm^2/(xm^2 +zm^2)^2))-(1/pi)*(3*(xm*(z_prime*zp + xm^2)/(xm^2 + zp^2)^2)-(3*(z_prime)^2*xm*zp*2 +xm^3*(4*z_prime^2 + 6*z*z_prime + z^2 + xm^2))/(xm^2+zp^2)^3);
Gxv= @(xm,zm,zp) (-1/(4*pi))*((zp/(xm^2 + zp^2))+ 2*((xm^2*zm/(xm^2+zm^2)^2)-(xm^2*zm)/(xm^2 +zm^2)^2))-(1/(2*pi)*(2*(zp/(xm^2 + zp^2)))-((2*z-z_prime)*(zp^2-xm^2)/(xm^2+zp^2)^2)+(2*z*z_prime*zp*(3*xm^2-zp^2))/(xm^2+zp^2)^3);
p= @(t) P_l*exp(-(-v*t)/2*alpha)/(4*pi*k*(-v*t));
term1 = -(alpha * E / (1 - 2*v)) * ...
    int(int((Gxh * dT_dx + Gxv * dT_dz), x_prime, -inf, inf), z_prime, 0, inf);
term2 = (2 * z) / pi * ...
    int(p(t) * (t - x)^2 / ((t - x)^2 + z^2)^2,t, -inf, inf);
term3 = -(alpha * E * T) / (1 - 2*v);
% Combine the terms
Sigma = term1 + term2 + term3;
% You can simplify Sigma if desired
simplifiedSigma = simplify(Sigma);
substitutedSigma=subs(simplifiedSigma,[t,x,z],[0,0.001,0]);

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Dyuman Joshi 2023-8-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2011692-matlab-not-computing-integral-of-an-infinite-integral#answer_1291957

在 MATLAB Online 中打开

Convert the symbolic functions to function handles, and use numerical integrals -

P_l=50;
v=0.1;
k=15;
Tm=1;
T0=300;
alpha=3.75*10^(-6);
E= 190*10^9;
nu=0.3;
syms x z x_prime z_prime t dT_dx dT_dz;
xm=x-x_prime;
zp=z+z_prime;
zm=z-z_prime;
% Define the terms
T=(P_l*exp(-((sqrt((x-v*t)^2 + (z)^2) + (x-v*t)))/(2*alpha)))/(4*3.14*k*sqrt((x-v*t)^2 + (z)^2)) +T0;
dT_dx=diff(T,x);
dT_dx_prime=subs(dT_dx,[x,z],[x_prime,z_prime]);
dT_dz=diff(T,z);
dT_dz_prime=subs(dT_dz,[x,z],[x_prime,z_prime]);
Gxh= (1/(4*pi))*(3*(xm/(xm^2 + zp^2)) + 2*(xm*zm^2/(xm^2 +zm^2)^2))-(1/pi)*(3*(xm*(z_prime*zp + xm^2)/(xm^2 + zp^2)^2)-(3*(z_prime)^2*xm*zp*2 +xm^3*(4*z_prime^2 + 6*z*z_prime + z^2 + xm^2))/(xm^2+zp^2)^3);
Gxv= (-1/(4*pi))*((zp/(xm^2 + zp^2))+ 2*((xm^2*zm/(xm^2+zm^2)^2)-(xm^2*zm)/(xm^2 +zm^2)^2))-(1/(2*pi)*(2*(zp/(xm^2 + zp^2)))-((2*z-z_prime)*(zp^2-xm^2)/(xm^2+zp^2)^2)+(2*z*z_prime*zp*(3*xm^2-zp^2))/(xm^2+zp^2)^3);
p= P_l*exp(-(-v*t)/2*alpha)/(4*pi*k*(-v*t));
%Convert to a function handle
T = matlabFunction(T);
p = matlabFunction(p);
fun = matlabFunction(Gxh * dT_dx + Gxv * dT_dz);
%Define terms as funciton handles
term1 = @(t,x,z) integral2(@(x_prime,z_prime) fun(t,x,z,x_prime,z_prime),-inf,inf,0,inf);
term2 = @(x,z) integral(@(t,x,z) (2.*z)./pi*(p(t).*(t - x).^2./((t - x).^2 + z^2).^2), t, -inf, inf);
term3 = @(t,x,z) -(alpha * E * T(t,x,z)) / (1 - 2*v);
Sigma = term1(0,0.001,0) + 0 + term3(0,0.001,0)
Sigma = -2.6719e+08

2 个评论
显示无隐藏无

AD 2023-8-23

Heyy..I am getting an error while calculating the term 2. ALso, there are no changes in the result for different t,x,z values.

Dyuman Joshi 2023-8-28

编辑：Dyuman Joshi 2023-8-28

在 MATLAB Online 中打开

I corrected the error for term2.

P_l=50;
v=0.1;
k=15;
Tm=1;
T0=300;
alpha=3.75*10^(-6);
E= 190*10^9;
nu=0.3;
syms x z x_prime z_prime t dT_dx dT_dz;
xm=x-x_prime;
zp=z+z_prime;
zm=z-z_prime;
% Define the terms
T=(P_l*exp(-((sqrt((x-v*t)^2 + (z)^2) + (x-v*t)))/(2*alpha)))/(4*3.14*k*sqrt((x-v*t)^2 + (z)^2)) +T0;
dT_dx=diff(T,x);
dT_dx_prime=subs(dT_dx,[x,z],[x_prime,z_prime]);
dT_dz=diff(T,z);
dT_dz_prime=subs(dT_dz,[x,z],[x_prime,z_prime]);
Gxh= (1/(4*pi))*(3*(xm/(xm^2 + zp^2)) + 2*(xm*zm^2/(xm^2 +zm^2)^2))-(1/pi)*(3*(xm*(z_prime*zp + xm^2)/(xm^2 + zp^2)^2)-(3*(z_prime)^2*xm*zp*2 +xm^3*(4*z_prime^2 + 6*z*z_prime + z^2 + xm^2))/(xm^2+zp^2)^3);
Gxv= (-1/(4*pi))*((zp/(xm^2 + zp^2))+ 2*((xm^2*zm/(xm^2+zm^2)^2)-(xm^2*zm)/(xm^2 +zm^2)^2))-(1/(2*pi)*(2*(zp/(xm^2 + zp^2)))-((2*z-z_prime)*(zp^2-xm^2)/(xm^2+zp^2)^2)+(2*z*z_prime*zp*(3*xm^2-zp^2))/(xm^2+zp^2)^3);
p= P_l*exp(-(-v*t)/2*alpha)/(4*pi*k*(-v*t));
%Convert to a function handle
T0 = matlabFunction(T);
p = matlabFunction(p);
f = Gxh * dT_dx + Gxv * dT_dz;
fun = matlabFunction(f);
%Define terms as function handles
term1 = @(t,x,z) integral2(@(x_prime,z_prime) fun(t,x,z,x_prime,z_prime), -inf, inf, 0, inf);
term2 = @(x,z) integral(@(t) (2.*z)./pi*(p(t).*(t - x).^2./((t - x).^2 + z^2).^2), -inf, inf);
term3 = @(t,x,z) -(alpha * E * T0(t,x,z)) / (1 - 2*v);

"ALso, there are no changes in the result for different t,x,z values."

Because the result is dominated by term3, in which there is not much change w.r.t values

format long
%t x z values
%0 0.001 0
term1(0,0.001,0)
ans = 
   -3.453931570108182e-107
term2(0.001,0)
ans = 
     0
term3(0,0.001,0)
ans = 
    -2.671874999999999e+08
%-0.5 0 5
term1(-0.5,0,5)
ans = 
     0
term2(0,5)
ans = 
    -4.973591972400049e-07
term3(-0.5,0,5)
ans = 
    -2.671874999999999e+08
%-5e3 0 0
term1(-5e3,0,0)
ans = 
     0
term2(0,0)
ans = 
     0
term3(-5e3,0,0)
ans = 
    -2.671874999999999e+08