trapezoidal rule for integral function

I need to use this rule for the following function, I tried but I obtained error, here the code and in the butom I show a function that I want to use there the trapezoidal

clear all;
xl=0; xr=1;            %domain[xl,xr]                                                                                                                         
J=100;                  % J: number of dividion for x
dx=(xr-xl)/J;          % dx: mesh size
tf=1000;                % final simulation time
Nt=100000;                 %Nt:number of time steps
dt=tf/Nt;                                                                                     
c=50;     
s=600;   % paremeter for the solution
lambdau1=0.0004;          %turning rate of susp(right)
lambdau2=0.0002;          %turning rate of susp(left)
lambdav1=0.0003;           %turning rate of infect(right)
lambdav2=0.0001;           %turning rate of infect(left)
beta=0.1;               % transimion rate
alpha=0.0;               %recovery/death rate
b=0.1;
d=0.5;
bV=0.3;
dV=0.1;
au=0.0003;                   % advection speed - susceptible
av=0.0001;                   % advection speed - infected
A=0.01;
mu1=au*dt/dx; 
mu2=av*dt/dx;
w = linspace(xl,xr,J);   
if mu1>1.0               %make sure dt satisfy stability condition
    error('mu1 should<1.0!')
end
if mu2>1.0               %make sure dt satisfy stability condition
    error('mu2 should<1.0!')
end
%Evaluate the initial conditions 
x=xl:dx:xr;               %generate the grid point
%fr=0*exp(-c*(x-0.5).^2);     %dimension f(1:J+1)initial conditions of right moving suscp
%fl=0*exp(-c*(x-0.5).^2);    %initial conditions of left moving suscep
%gr=0.5*exp(-s*(x-0.5).^2);  %initial conditions of infected
%gl=0.5*exp(-s*(x-0.5).^2);  %initial conditions of infected
k=6*pi/(xr-xl);
fr=0+0.01*(1+sin(k*x));     %dimension f(1:J+1)initial conditions of right moving suscp
fl=0+0.01*(1+sin(k*x));    %initial conditions of left moving suscep
gr=0.25+0.01*(1+sin(k*x));  %initial conditions of infected
gl=0.75+0.01*(1+sin(k*x));
%store the solution at all grid points for all time steps
ur=zeros(J+1,Nt);
ul=zeros(J+1,Nt);
vr=zeros(J+1,Nt);
vl=zeros(J+1,Nt);
v=zeros(J+1,Nt);
K=1/(0.05*sqrt(2*pi))*exp(-0.5*(s/0.05)^2);
y= b*K*v(x+s)+d*K*ur(x+s);

Array indices must be positive integers or logical values.

F = int(y,w);
z = trapz(w, y);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tt=0:tf/(Nt-1):tf;
[T,X]=meshgrid(tt,x);
%find the approximate solution at each time step 
for  n=1:Nt
    t=n*dt;               % current time
    if n==1                         % first time step
        for j=2:J                      % interior nodes
            %%%%%%%%%%%%%%%%%%%upwind of right-moving of susceptible
            ur(j,n)=fr(j)-0.5*mu1*(fr(j)-fr(j-1))+dt*lambdau2*(0.5+0.5*tanh(1-z))*fl(j)-dt*lambdau1*(0.5+0.5*tanh(1-z))*fr(j)-dt*beta*fr(j)*(gr(j)+gl(j))+dt*0.5*alpha*(gr(j)+gl(j));
            %%%%%%%%%%%%%%%%%%%upwind of left-moving of susceptible
            ul(j,n)=fl(j)+0.5*mu1*(fl(j+1)-fl(j))+dt*lambdau1*(0.5+0.5*tanh(1-z))*fr(j)-dt*lambdau2*(0.5+0.5*tanh(1-z))*fl(j)-dt*beta*fl(j)*(gr(j)+gl(j))+dt*0.5*alpha*(gr(j)+gl(j));
            %%%%%%%%%%%%%%%%%%%upwind of right-moving of infected
            vr(j,n)=gr(j)-0.5*mu2*(gr(j)-gr(j-1))+dt*lambdav2*(0.5+0.5*tanh(1-bV*(gr(j)+gl(j))+dV*fr(j)))*gl(j)-dt*lambdav1*(0.5+0.5*tanh(1-bV*(gr(j)+gl(j))+dV*fl(j)))*gr(j)+dt*beta*fr(j)*(gr(j)+gl(j))-dt*alpha*gr(j);
            %%%%%%%%%%%%%%%%%%%lupwind of left-moving of infected
            vl(j,n)=gl(j)+0.5*mu2*(gl(j+1)-gl(j))+dt*lambdav1*(0.5+0.5*tanh(1-bV*(gr(j)+gl(j))+dV*fl(j)))*gr(j)-dt*lambdav2*(0.5*+0.5*tanh(1-bV*(gr(j)+gl(j))+dV*fr(j)))*gl(j)+dt*beta*fl(j)*(gr(j)+gl(j))-dt*alpha*gl(j);
            
            
        end
        %peridic BC for right moving of suscep(j=J+1,j+2=1)     
        ur(J+1,1)= fr(J+1)-0.5*mu1*(fr(J+1)-fr(J))+dt*lambdau2*(0.5+0.5*tanh(1-b*(gr(J+1)+gl(J+1))+d*fr(J+1)))*fl(J+1)-dt*lambdau1*(0.5+0.5*tanh(1-b*(gr(J+1)+gl(J+1))+d*fl(J+1)))*fr(J+1)-dt*beta*fr(J+1)*(gr(J+1)+gl(J+1))+dt*0.5*alpha*(gr(J+1)+gl(J+1)); 
        %peridic BC for right moving of suscep(j=1,0=J)
        ur(1,1)=fr(1)-0.5*mu1*(fr(1)-fr(J))+dt*lambdau2*(0.5+0.5*tanh(1-b*(gr(1)+gl(1))+d*fr(1)))*fl(1)-dt*lambdau1*(0.5+0.5*tanh(1-b*(gr(1)+gl(1))+d*fl(1)))*fr(1)-dt*beta*fr(1)*(gr(1)+gl(1))+dt*0.5*alpha*(gr(1)+gl(1)); %peridic BC for right moving
        %peridic BC for left moving of suscep(j=1)
        ul(1,1)=fl(1)+0.5*mu1*(fl(2)-fl(1))+dt*lambdau1*(0.5+0.5*tanh(1-b*(gr(1)+gl(1))+d*fl(1)))*fr(1)-dt*lambdau2*(0.5*+0.5*tanh(1-b*(gr(1)+gl(1))+d*fr(1)))*fl(1)-dt*beta*fl(1)*(gr(1)+gl(1))+dt*0.5*alpha*(gr(1)+gl(1)); %peridic BC for left moving;  %peridic BC for left moving
        %peridic BC for left moving of suscep(j=J+1)
        ul(J+1,1)=fl(J+1)+0.5*mu1*(fl(1)-fl(J+1))+dt*lambdau1*(0.5+0.5*tanh(1-b*(gr(J+1)+gl(J+1))+d*fl(J+1)))*fr(J+1)-dt*lambdau2*(0.5*+0.5*tanh(1-b*(gr(J+1)+gl(J+1))+d*fr(J+1)))*fl(J+1)-dt*beta*fl(J+1)*(gr(J+1)+gl(J+1))+dt*0.5*alpha*(gr(J+1)+gl(J+1));  %peridic BC for right moving(j=J+1);  %peridic BC for left moving
        %peridic BC for right moving of infected(j=J+1)
        vr(J+1,n)=gr(J+1)-0.5*mu2*(gr(J+1)-gr(J))+dt*lambdav2*(0.5+0.5*tanh(1-bV*(gr(J+1)+gl(J+1))+dV*fr(J+1)))*gl(J+1)-dt*lambdav1*(0.5+0.5*tanh(1-bV*(gr(J+1)+gl(J+1))+dV*fl(J+1)))*gr(J+1)+dt*beta*fr(J+1)*(gr(J+1)+gl(J+1))-dt*alpha*gr(J+1);
        %peridic BC for right moving of infected(j=1)
        vr(1,n)=gr(1)-0.5*mu2*(gr(1)-gr(J))+dt*lambdav2*(0.5+0.5*tanh(1-bV*(gr(1)+gl(1))+dV*fr(1)))*gl(1)-dt*lambdav1*(0.5+0.5*tanh(1-bV*(gr(1)+gl(1))+dV*fl(1)))*gr(1)+dt*beta*gr(1)*(gr(1)+gl(1))-dt*alpha*gr(1);
        %peridic BC for left moving of infected(j=1)
        vl(1,n)=gl(1)+0.5*mu2*(gl(2)-gl(1))+dt*lambdav1*(0.5+0.5*tanh(1-bV*(gr(1)+gl(1))+dV*fl(1)))*gr(1)-dt*lambdav2*(0.5*+0.5*tanh(1-bV*(gr(1)+gl(1))+dV*fr(1)))*gl(1)+dt*beta*fl(J+1)*(gr(1)+gl(1))-dt*alpha*gl(1);
        %peridic BC for left moving of infected(j=J+1)
        vl(J+1,n)=gl(J+1)+0.5*mu2*(gl(1)-gl(J+1))+dt*lambdav1*(0.5+0.5*tanh(1-bV*(gr(J+1)+gl(J+1))+dV*fl(J+1)))*gr(J+1)-dt*lambdav2*(0.5*+0.5*tanh(1-bV*(gr(J+1)+gl(J+1))+dV*fr(J+1)))*gl(J+1)+dt*beta*fl(J+1)*(gr(J+1)+gl(J+1))-dt*alpha*gl(J+1);
    else
        for j=2:J                % interior nodes
            %%%%%%%%%%%%%%%%%%%upwind of right-moving of susceptible
            ur(j,n)=ur(j,n-1)-0.5*mu1*(ur(j,n-1)-ur(j-1,n-1))+dt*lambdau2*(0.5+0.5*tanh(1-b*(vr(j,n-1)+vl(j,n-1))+d*ur(j,n-1)))*ul(j,n-1)-dt*lambdau1*(0.5+0.5*tanh(1-b*(vr(j,n-1)+vl(j,n-1))+d*ul(j,n-1)))*ur(j,n-1)-dt*beta*ur(j,n-1)*(vr(j,n-1)+vl(j,n-1))+dt*0.5*alpha*(vr(j,n-1)+vl(j,n-1));
            %%%%%%%%%%%%%%%%%%%upwind of right-moving of susceptible
            ul(j,n)=ul(j,n-1)+0.5*mu1*(ul(j+1,n-1)-ul(j,n-1))+dt*lambdau1*(0.5+0.5*tanh(1-b*(vr(j,n-1)+vl(j,n-1))+d*ul(j,n-1)))*ur(j,n-1)-dt*lambdau2*(0.5*+0.5*tanh(1-b*(vr(j,n-1)+vl(j,n-1))+d*ur(j,n-1)))*ul(j,n-1)-dt*beta*ul(j,n-1)*(vr(j,n-1)+vl(j,n-1))+dt*0.5*alpha*(vr(j,n-1)+vl(j,n-1));    
            %%%%%%%%%%%%%%%%%%%upwindof right-moving of infected
            vr(j,n)=vr(j,n-1)-0.5*mu2*(vr(j,n-1)-vr(j-1,n-1))+dt*lambdav2*(0.5+0.5*tanh(1-bV*(vr(j,n-1)+vl(j,n-1))+dV*ur(j,n-1)))*vl(j,n-1)-dt*lambdav1*(0.5+0.5*tanh(1-bV*(vr(j,n-1)+vl(j,n-1))+dV*ul(j,n-1)))*vr(j,n-1)+dt*beta*ur(j,n-1)*(vr(j,n-1)+vl(j,n-1))-dt*alpha*vr(j,n-1);
            %%%%%%%%%%%%%%%%%%%upwind of left-moving of infected
            vl(j,n)=vl(j,n-1)+0.5*mu2*(vl(j+1,n-1)-vl(j,n-1))+dt*lambdav1*(0.5+0.5*tanh(1-bV*(vr(j,n-1)+vl(j,n-1))+dV*ul(j,n-1)))*vr(j,n-1)-dt*lambdav2*(0.5*+0.5*tanh(1-bV*(vr(j,n-1)+vl(j,n-1))+dV*ur(j,n-1)))*vl(j,n-1)+dt*beta*ul(j,n-1)*(vr(j,n-1)+vl(j,n-1))-dt*alpha*vl(j,n-1);
            
        end 
        %peridic BC for right moving of suscep(j=J+1,j+2=1)
        ur(J+1,n)=ur(J+1,n-1)-0.5*mu1*(ur(J+1,n-1)-ur(J,n-1))+dt*lambdau2*(0.5+0.5*tanh(1-b*(vr(J+1,n-1)+vl(J+1,n-1))+d*ur(J+1,n-1)))*ul(J+1,n-1)-dt*lambdau1*(0.5+0.5*tanh(1-b*(vr(J+1,n-1)+vl(J+1,n-1))+d*ul(J+1,n-1)))*ur(J+1,n-1)-dt*beta*ur(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))+dt*0.5*alpha*(vr(J+1,n-1)+vl(J+1,n-1));
        %peridic BC for right moving of suscep(j=1,0=J)
        ur(1,n)=ur(1,n-1)-0.5*mu1*(ur(1,n-1)-ur(J,n-1))+dt*lambdau2*(0.5+0.5*tanh(1-b*(vr(1,n-1)+vl(1,n-1))+d*ur(1,n-1)))*ul(1,n-1)-dt*lambdau1*(0.5+0.5*tanh(1-b*(vr(1,n-1)+vl(1,n-1))+d*ul(1,n-1)))*ur(1,n-1)-dt*beta*ur(1,n-1)*(vr(1,n-1)+vl(1,n-1))+dt*0.5*alpha*(vr(1,n-1)+vl(1,n-1));
        %peridic BC for left moving of suscep (j=1)
        ul(1,n)=ul(1,n-1)+0.5*mu1*(ul(2,n-1)-ul(1,n-1))+dt*lambdau1*(0.5+0.5*tanh(1-b*(vr(1,n-1)+vl(1,n-1))+d*ul(1,n-1)))*ur(1,n-1)-dt*lambdau2*(0.5*+0.5*tanh(1-b*(vr(1,n-1)+vl(1,n-1))+d*ur(1,n-1)))*ul(1,n-1)-dt*beta*ul(1,n-1)*(vr(1,n-1)+vl(1,n-1))+dt*0.5*alpha*(vr(1,n-1)+vl(1,n-1));    
        %peridic BC for left moving of suscep(j=J+1)
        ul(J+1,n)=ul(J+1,n-1)+0.5*mu1*(ul(1,n-1)-ul(J+1,n-1))+dt*lambdau1*(0.5+0.5*tanh(1-b*(vr(J+1,n-1)+vl(J+1,n-1))+d*ul(J+1,n-1)))*ur(J+1,n-1)-dt*lambdau2*(0.5*+0.5*tanh(1-b*(vr(J+1,n-1)+vl(J+1,n-1))+d*ur(J+1,n-1)))*ul(J+1,n-1)-dt*beta*ul(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))+dt*0.5*alpha*(vr(J+1,n-1)+vl(J+1,n-1));    
        %peridic BC for right moving of infected(j=J+1)
        vr(J+1,n)=vr(J+1,n-1)-0.5*mu2*(vr(J+1,n-1)-vr(J,n-1))+dt*lambdav2*(0.5+0.5*tanh(1-bV*(vr(J+1,n-1)+vl(J+1,n-1))+dV*ur(J+1,n-1)))*vl(J+1,n-1)-dt*lambdav1*(0.5+0.5*tanh(1-bV*(vr(J+1,n-1)+vl(J+1,n-1))+dV*ul(J+1,n-1)))*vr(J+1,n-1)+dt*beta*ur(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))-dt*alpha*vr(J+1,n-1);
        %peridic BC for right moving of infected(j=1)
        vr(1,n)=vr(1,n-1)-0.5*mu2*(vr(1,n-1)-vr(J,n-1))+dt*lambdav2*(0.5+0.5*tanh(1-bV*(vr(1,n-1)+vl(1,n-1))+dV*ur(1,n-1)))*vl(1,n-1)-dt*lambdav1*(0.5+0.5*tanh(1-bV*(vr(1,n-1)+vl(1,n-1))+dV*ul(1,n-1)))*vr(1,n-1)+dt*beta*ur(1,n-1)*(vr(1,n-1)+vl(1,n-1))-dt*alpha*vr(1,n-1);
        %peridic BC for left moving of infected(j=1)
        vl(1,n)=vl(1,n-1)+0.5*mu2*(vl(2,n-1)-vl(1,n-1))+dt*lambdav1*(0.5+0.5*tanh(1-bV*(vr(1,n-1)+vl(1,n-1))+dV*ul(1,n-1)))*vr(1,n-1)-dt*lambdav2*(0.5*+0.5*tanh(1-bV*(vr(1,n-1)+vl(1,n-1))+dV*ur(1,n-1)))*vl(1,n-1)+dt*beta*ul(1,n-1)*(vr(1,n-1)+vl(1,n-1))-dt*alpha*vl(1,n-1);
        %peridic BC for left moving of infected(j=J+1)
        vl(J+1,n)=vl(J+1,n-1)+0.5*mu2*(vl(1,n-1)-vl(J+1,n-1))+dt*lambdav1*(0.5+0.5*tanh(1-bV*(vr(J+1,n-1)+vl(J+1,n-1))+dV*ul(J+1,n-1)))*vr(J+1,n-1)-dt*lambdav2*(0.5*+0.5*tanh(1-bV*(vr(J+1,n-1)+vl(J+1,n-1))+dV*ur(J+1,n-1)))*vl(J+1,n-1)+dt*beta*ul(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))-dt*alpha*vl(J+1,n-1);
    end
    
    %calculate the analytical solution 
    for j=1:J+1
        xj=xl+(j-1)*dx;
        u_ex(j,n)=exp(-c*(xj-au*t-0.5)^2);
    end
end 
figure
ss=surf(T,X,ur+ul);
set(ss,'edgecolor','none');
xlabel('t','FontSize',28);
ylabel('x','FontSize',28);
title('u^++u^-','FontSize',28);
set(gca,'Fontsize',28);
view(2)
%%%%%%%%%%%%
figure
ss=surf(T,X,vr+vl);
set(ss,'edgecolor','none');
xlabel('t','FontSize',28);
ylabel('x','FontSize',28);
title('v^++v^-','FontSize',28);
set(gca,'Fontsize',28);
view(2)
figure
ss=surf(T,X,ur);
set(ss,'edgecolor','none');
xlabel('t','FontSize',28);
ylabel('x','FontSize',28);
title('u^+','FontSize',28);
set(gca,'Fontsize',28);
figure
ss=surf(T,X,vr);
set(ss,'edgecolor','none');
xlabel('t','FontSize',28);
ylabel('x','FontSize',28);
title('v^+','FontSize',28);
set(gca,'Fontsize',28);
figure
ss=surf(T,X,ul);
set(ss,'edgecolor','none');
xlabel('t','FontSize',28);
ylabel('x','FontSize',28);
title('u^-','FontSize',28);
set(gca,'Fontsize',28);
figure
ss=surf(T,X,vl);
set(ss,'edgecolor','none');
xlabel('t','FontSize',28);
ylabel('x','FontSize',28);
title('v^-','FontSize',28);
set(gca,'Fontsize',28);
figure
ss=surf(T,X,vr+vl);
set(ss,'edgecolor','none');
xlabel('t','FontSize',28);
ylabel('x','FontSize',28);
title('v^++v^-','FontSize',28);
set(gca,'Fontsize',28);
figure
ss=surf(T,X,ur+ul);
set(ss,'edgecolor','none');
xlabel('t','FontSize',28);
ylabel('x','FontSize',28);
title('u^++u^-','FontSize',28);
set(gca,'Fontsize',28); 

Answer 1

[Akshat Dalal]

Hi lulu,

I understand you are facing some error while implementing the trapezoidal rule for integration. The error seems to be happening in the following line of code:

y= b*K*v(x+s)+d*K*ur(x+s);

Here, the variable ‘x’ is an array with non-integral values, while 's' is a constant integer. Thus, the array ‘x + s’ will contain some non-integral values which cannot be used as array indices and throws an error.