# MATLAB integral3 Error

1 view (last 30 days)
Shahid Hasnain on 19 Nov 2018
Here is my Code, dexact_v3 is not working with some issue. Your help will highly appreciated.
t=0;L = 1; hx=0.01;
a=0;b=L;
x=a:hx:b;y=x;
My_result_W=wexact_Patrick(t,x,y);
%My_result_d=dexact_v3(t,x,y);
function [d1]=dexact_v3(t,x,y);
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);
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);
% The final result for d1 is
return
%%%%%%%%%%%%%%%%%%%%%%% Working Code which I used in above Code%%%%%%%%%
function [W1]=wexact_Patrick(t,x,y);
L=1;
Ntrunc = 10;
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
for p = 1:Ntrunc
for q = 1: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);
end
end
% The analytical solution for w(x,y,t) and first order of d(x,y,t)
W1 = 1; %Constant term
for p =1: Ntrunc
for q =1: Ntrunc
r = p^2 + q^2 + 1;
W1 = W1 + bd(p,q).*exp(-r.*t).*sin(p.*x).*cos(q.*y);
end
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%