Long Numerical Integration: Is this possible with MATLAB?

Join the discussion
Ask a follow-up here, or get your own question answered by working scientists, mathematicians and engineers — people, not an autocomplete.
Real named experts · corrections over time · the nuance an AI answer skips
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.
 
Thanks, let me try it.