trapezoidal rule for integral function

2 次查看（过去 30 天）
lulu
lulu 2022-8-22
编辑： Akshat Dalal 2023-10-8
Ran in:
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);

请先登录，再回答此问题。

回答（1 个）

Akshat Dalal
Akshat Dalal 2023-10-4
编辑：Akshat Dalal 2023-10-8
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.

类别

在 Help Center 和 File Exchange 中查找有关 Spline Postprocessing 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by