Numerical integration of pressure response with infinite sum in MATLAB

Join the discussion
Registration is free. Start your own thread to ask a follow-up.
4 replies · 3K views
Edgarg
Messages
3
Reaction score
0
Please, I am new in MATLAB and need some help as to whether MATLAB can perform this integration, and how to go about it. I have tried 'quad', 'quad8', 'trapz' and even a sum approach but it returns either a not-a-number answer or mtimes error or mrdivide error message. What I want to do is this (Pressure response of a horizontal observation well at a distance xd, yd and zd from a producer)
integral= (1/(sqrt(td)))*(exp(yd^2/(4*td)))*(erf((1+xd)/(2*sqrt(td)))+erf((1-xd)/(2*sqrt(td))))*(1+2*(exp((-n^2*pi^2*td)/hd^2)*cos(n*pi*zwd)*cos(n*pi*((zd/hd)+zwd))));
xd, yd, zd, hd and zwd are constants.
i am integrating from 0 to td for td values of 0.001 to 10000
for each of the integral segment, 0.001 to 10000, (exp((-n^2*pi^2*td)/hd^2)*cos(n*pi*zwd)*cos(n*pi*((zd/hd)+zwd))) will be summed for n=1 to infinity. ie Ʃ(exp((-n^2*pi^2*td)/hd^2)*cos(n*pi*zwd)*cos(n*pi*((zd/hd)+zwd))) for n=1 to ∞

This is my recent attempt: where I did td =1 to 10000; and n= 1 to 100 just for it to work first.
----------------------
clc
x=2800;
y=0;
z=50;
zw=50.165;
kx=300;
ky=300;
kz=100;
l=1000; %well half length
h=100;
% t=150;
% td=(0.0002637*ky*t)/(phi*mu*ct*xf.^2)
k=nthroot(kx*ky*kz,3);
xd = x/l*sqrt(k/kx);
yd = y/l*sqrt(k/ky);
zd = z/h;
zwd=zw/h;
hd = h/l*sqrt(k/kz);
for i = 1:1e+1:10000
for j=1:100
int(i,j)= (1./(sqrt(td(i)))).*(exp(yd.^2./(4.*td(i)))).*(erf((1+xd)./(2.*sqrt(td(i))))+erf((1-xd)./(2.*sqrt(td(i))))).*(1+2.*(exp((-n(j)^2.*pi^2.*td(i))./hd.^2).*cos(n(j)*pi*zwd).*cos(n(j)*pi*((zd./hd)+zwd))));
end
end
p=sum(int*10,3)

Error: ? Undefined function or method 'td' for input arguments of type 'double'.
**The attached document contains the equation.
 
Physics news on Phys.org
I'm having trouble understanding this. Is this all of your code? You don't have a definition for td.
 
Yes, I am integrating with respect to td. So, I did not define td.
Lower integral limit is zero, and I will have several upper limits of 0.001 to 10000 with varying intervals. The equation contains a sum series within.
Equation is:
integral= (1/(sqrt(td)))*(exp(yd^2/(4*td)))*(erf((1+xd)/(2*sqrt(td)))+erf((1-xd)/(2*sqrt(td))))*(1+2*Ʃ(exp((-n^2*pi^2*td)/hd^2)*cos(n*pi*zwd)*cos(n*pi*((zd/hd)+zwd))))
the sigma is for n=1 to infinity.
 
If you're using it as a symbolic variable, you need to define it using
Code:
syms td;

I don't know if there are alternate syntaxes for integrating, but AFAIK it's normally

Code:
int(expr,var)

so your code would be

Code:
syms td
int((1./(sqrt(td(i)))).*(exp(yd.^2./(4.*td(i)))).*(erf((1+xd)./(2.*sqrt(td(i))))+erf((1-xd)./(2.*sqrt(td(i))))).*(1+2.*(exp((-n(j)^2.*pi^2.*td(i))./hd.^2).*cos(n(j)*pi*zwd).*cos(n(j)*pi*((zd./hd)+zwd)))),td)

or something like that. If you're trying to integrate numerically rather than symbolically, then you should evaluate the function first and then integrate it, and you would need to define td.