Error using integral2

I have following expression with

x=[1 2 3 4 5]
y=[1 2 3 4 5]
L = 100; % Length of the box
p=1;
q=1;
Ui = @(x,y)  sin(pi*x).*sin(pi*y); % An initial distribution for u, i.e. the u(x,y,0)
Vi = @(x,y)  sin(pi*x).*sin(pi*y); % An initial distribution for v, i.e. the v(x,y,0)
Wi = Ui(x,y) + Vi(x,y); % Initial condition for w, i.e. w(x,y,0)
% Computation of the b-coefficients
if p < Ntrunc
	if q < Ntrunc
		my_integrand = @(x,y)  4*Wi.*sin(p.*x).*cos(q.*y)/L^2; 
        bd(p,q) = integral2(my_integrand(x,y),0,L,0,L);
		q = q+1;
	end
	p = p+1;
end 

I got following error,

Error using integral2Calc>integral2t/tensor (line 231)

Input function must return 'double' or 'single' values. Found 'function_handle'.

Error in integral2Calc>integral2t (line 55)

[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in integral2Calc (line 9)

[q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral2 (line 106)

Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);

Correction will be highly appreciated.

[madhan ravi]

Ntrunc?

[Shahid Hasnain]

Ntruc=100;

number of trunction terms. We say as 100.

[madhan ravi]

see my answer below

Answer 1

[Star Strider]

This works, although it throws a warning:

Ntrunc = 100;
L = 100; % Length of the box
p=1;
q=1;
Ui = @(x,y)  sin(pi*x).*sin(pi*y); % An initial distribution for u, i.e. the u(x,y,0)
Vi = @(x,y)  sin(pi*x).*sin(pi*y); % An initial distribution for v, i.e. the v(x,y,0)
Wi = @(x,y) Ui(x,y) + Vi(x,y); % Initial condition for w, i.e. w(x,y,0)
% Computation of the b-coefficients
if p < Ntrunc
	if q < Ntrunc
		my_integrand = @(x,y)  4*Wi(x,y).*sin(p.*x).*cos(q.*y)/L^2; 
        bd(p,q) = integral2(@(x,y)my_integrand(x,y),0,L,0,L);
		q = q+1;
	end
	p = p+1;
end 

I defer to you to correct the problem that produces the warning.

EDIT — Corrected typographical error.

[Shahid Hasnain]

@Star Strider

Thank you and highly appreciated your answer. It is working now. Let me check other files also.

Again Thank you.

[Shahid Hasnain]

Again half code is working and half got error .Let me post complete code here as

t=0; x=[1 2 3 4 5]; y=[1 2 3 4 5];

function [W1]=wexact_Patrick(t,x,y);

Ntrunc = 100;

L = 100; % Length of the box

p=1;

q=1;

Ui = @(x,y) sin(pi*x).*sin(pi*y); % An initial distribution for u, i.e. the u(x,y,0)

Vi = @(x,y) sin(pi*x).*sin(pi*y); % An initial distribution for v, i.e. the v(x,y,0)

Wi = @(x,y) Ui(x,y) + Vi(x,y); % Initial condition for w, i.e. w(x,y,0)

% Computation of the b-coefficients

if p < Ntrunc

if q < Ntrunc

my_integrand = @(x,y) 4*Wi(x,y).*sin(p.*x).*cos(q.*y)/L^2;

bd(p,q) = integral2(@(x,y)my_integrand(x,y),0,L,0,L);

q = q+1;

end

p = p+1;

end

% The analytical solution for w(x,y,t) and first order of d(x,y,t)

W1 = 1 %Constant term

if p < Ntrunc

if q < Ntrunc

r = p^2 + q^2 + 1;

W1 = W1 + bd(p,q).*exp(-r.*t).*sin(p.*x).*cos(q.*y);

q = q+1;

end

p = p+1;

end

return

Error is about

wexact_Patrick(t,x,y)

Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test.

> In integral2Calc>integral2t (line 129)

In integral2Calc (line 9)

In integral2 (line 106)

In wexact_Patrick (line 13)

W1 =

1

Attempted to access bd(2,2); index out of bounds because numel(bd)=1.

Error in wexact_Patrick (line 24)

W1 = W1 + bd(p,q).*exp(-r.*t).*sin(p.*x).*cos(q.*y);

Thank you in advance for correction.

[Star Strider]

I am not certain what you are doing, so I do not have any idea of how to suppress the Warning. My impression is that integral2 is having a difficult time with the ‘my_integrand’ function you are giving it. It may be encountering a discontinuity, although I did not plot it to see what it is doing. I must leave that to you to correct.

I changed the nested while loops into for loops, and made the evaluation marginally more efficient by moving ‘my_integrand’ out of the loop. My changes are specifically:

my_integrand = @(x,y,p,q)  4*Wi(x,y).*sin(p.*x).*cos(q.*y)/L^2; 
% Computation of the b-coefficients
for p = 1 : Ntrunc
	for q = 1 : Ntrunc
        bd(p,q) = integral2(@(x,y)my_integrand(x,y,p,q),0,L,0,L);
    end
end 

This creates the appropriate matrix, so it should eliminate this error:

Attempted to access bd(2,2); index out of bounds because numel(bd)=1.
Error in wexact_Patrick (line 24)
		W1 = W1 + bd(p,q).*exp(-r.*t).*sin(p.*x).*cos(q.*y);

I apologise for the delay. It was overnight for me here, and additionally, I could not figure out how to eliminate the Warning. I suspect it is due to the increasing magnitudes of ‘p’ and ‘q’, since it does not appear until they get reasonably large, and setting ‘Ntrunc’ to a relatively low value prevents the Warning altogether.

[Shahid Hasnain]

@Star Strider

Thank you, you did great, I fixed it with your great support.

Now we have another file. I do my best to explain as we have 3 files

1- [W1]=wexact_Patrick(t,x,y); which is now working. and I attached.

2-[d1]=dexact_v3(t,x,y); partially working upto line 36 and rest of the code I do comment which is also attached.

3- To test these I used one test file which is also attached. (comment on 10th line)

Now Problem in dexact_v3(t,x,y); file at the end. Please is it possible to rectify error?

Say highly thanks and great appreciations.

Thank you in advance again.

[Star Strider]

Please describe the problems you are having with ‘dexact_v3’.

Also, first please re-write your code to eliminate all the global variables. Then test it again with the global calls removed. You need to include those variables as extra parameters to your functions, instead. It is best to completely avoid global variables in all code.

[Shahid Hasnain]

@Star Strider,

Sorry for being late, Now this is working until line 42 (I paste it here)

%%%%%%%%%%%

t=0;

L = 1; % Length of the box

hx=0.01;

a=0;

b=L;

x=a:hx:b;

y=x;

Ntrunc = 10; % Number of summations (high frequency oscillations can be neglected)

T = 10; % Duration of the time interval

p=1;

q=1;

Ui = @(x,y) (cos(pi.*x).*cos(pi.*y)); % An initial distribution for u, i.e. the u(x,y,0)

Vi = @(x,y) (sin(pi.*x).*sin(pi.*y)); % An initial distribution for v, i.e. the v(x,y,0)

di = @(x,y) Ui(x,y) - Vi(x,y); % Initial condition for w, i.e. w(x,y,0)

% Computation of the b'-coefficients (denoted by bd(p,q))

for p=1:Ntrunc

for q =1:Ntrunc

u_integrand = @(x,y) 4*di(x,y).*sin(p.*x).*cos(q.*y)/L.^2;

m_bd(p,q)= integral2(@(x,y)u_integrand(x,y),0,L,0,L);

end

end

% Zeroth order of d(x,y,t) without the convolution integral term, this is denoted by dzeroh

dzeroh = 1; % Constant term

for p=1: Ntrunc

for q= 1:Ntrunc

r = p.^2 + q.^2 + 1;

dzeroh = @(x,y,t) (dzeroh(x,y,t) + m_bd(p,q).*exp(-r.*t).*sin(p.*x).*cos(q.*y));

q = q+1;

end

p = p+1;

end

hx = 0.1; % discretization of the x-interval

hy = 0.1; % " " " y- "

ht = 0.1; % " " " z- "

xc = a:hx:L;

yc = xc;

tc = 0:ht:T;

G= @(x,y,t,xc,yc,tc) ( sqrt(pi/(t-tc)).*e^(-(t-tc)).*exp(-((x-xc).^2+(y-yc).^2)/(4*(t-tc))) ); % The Green function

integrand = @(x,y,t,xc,yc,tc) (-1.*((wexact_Patrick(tc,xc,tc)+1).*(wexact_Patrick(tc,xc,yc)-1).^2).*G(x,y,t,xc,yc,tc)/4);

%%%%%%%%%% Working upto here%%%%%%%%%%%%%%

[Shahid Hasnain]

If I add the following expressions to above Code, I face problem,

dzeroadd = integral3(@(xc,yc,tc)integrand(xc,yc,tc),0,L,0,L,0,T); % Convolution term

% % % Add the convolution integral term to the homogeneous solution (the solution without the nonlinear term)

dzero = dzeroh + dzeroadd; % It holds: dzero = 1 + d'_0

% % Here the complete solution for d is computed

integrand = @(x,y,t,xc,yc,tc) (-1.*(wexact_Patrick(xc,yc,tc)+dzero(xc,yc,tc)).*(wexact_Patrick(xc,yc,tc)-dzero(xc,yc,tc)).^2.*G(x,y,t,xc,yc,tc)/4);

dadd = integral3(@(xc,yc,tc) integrand,0,L,0,L,0,T); % Another convolution term to be added

% The final result for d1 is

d1 = dzeroh + dadd;

[Star Strider]

I do not understand.

How does this all fit together? What problem are you having?

[Shahid Hasnain]

@Star Strider,

Sorry for being late.

I am trying to solve Gray Scott Model in 2 dimension by analytical technique. the code I attached few days back, it was a try to find analytical solution.

Here is the system

u_t = u_{xx} + u_{yy} - uv^2 + 1 - u

v_t = v_{xx} + v_{yy} + uv^2 - v

I do not have analytical solution to this problem, thats why we choose another way to find.

w=u+v (wexact_Patrick, it was mentioned in previous posts) working

d=u-v (it was mentioned in previous posts) half working still.

Still working to rearrange integral3 which create problems.

[Star Strider]

No worries about a delayed reply.

Do you want to integrate a system of partial differential equations? If so, please post this as a new Question. My experience with the Partial Diufferential Equation Toolbox is limited, although Torsten and some others here likely have significant experience with it. I have not worked with it in a while, so it would take some time for me to provide you with a response.