Clear Filters
Clear Filters

MATLAB integral3 Error

2 views (last 30 days)
Shahid Hasnain
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);
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;
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Answers (0)

Categories

Find more on MATLAB in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!