How to do numerical integration using trapezoidal rule?

7 visualizzazioni (ultimi 30 giorni)
I want a graph for velocity profile. In my equation have one of the term is integro_differential. for that purpose i have to solve the integral part in my equation using newtons cotes closed integration. so, I'm using one of the rule is trapezoidal form. I want to use trapezoidal rule for integration. i.e . Here, we have to notice the 'T' is another equation output value which I already solved and i got the matrix form. i.e
%===========================================%
ymax=20; m=80; dy=ymax/m; %y=dy:dy:ymax;
xmax=1; n=20; dx=xmax/n;
tmax=100; nt=500; dt=tmax/nt; t=0:dt:tmax;
tol=1e-2;
max_difference(1)=1; max_Iteration=1;
pr=6.2;
UOLD=zeros(m,nt);VOLD=zeros(m,nt);
TNEW=0; TOLD=TNEW*ones(m,nt); TWALL=ones(1,length(t));
A=zeros([1,m]);
B=A;
C=A;
D=A;
T=TOLD;
while max_Iteration>tol
for j=1:nt
for i=1:m
if j>i
C(i)=dt*VOLD(i,j)/4*dy-dt/(2*dy^2);
elseif i>j
A(i)=-dt*VOLD(i,j)/4*dy-dt/(2*dy^2);
elseif i==j
B(i)=1+dt*UOLD(i,j)/2*dx+dt/dy^2;
end
end
end
for j=2:nt
for i=2:m-1
if i==2
D(i)=-(dt*UOLD(i,j)*(-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1)/2*dx))+(dt*(TWALL(j)-2*TOLD(i,j)+TOLD(i-1,j)/(2*dy^2)))-(dt*VOLD(i,j)*(TOLD(i-1,j)-TWALL(j)+TOLD(i,j)/4*dy))-((-dt)*(VOLD(i,j)/4*dy))-(dt*(TWALL(j)/2*dy^2));
elseif i==m-1
D(i)=-(dt*UOLD(i,j)*(-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1)/2*dx))+(dt*(TOLD(i+1,j)-2*TOLD(i,j)+TNEW/(2*dy^2)))-(dt*VOLD(i,j)*(TNEW-TOLD(i+1,j)+TOLD(i,j)/4*dy))-(dt*(VOLD(i,j)/4*dy))-(dt*(TNEW/2*dy^2));
else
D(i)=-(dt*UOLD(i,j)*(-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1)/2*dx))+(dt*(TOLD(i+1,j)-2*TOLD(i,j)+TOLD(i-1,j)/2*dy^2))-(dt*VOLD(i,j)*(TOLD(i-1,j)-TOLD(i+1,j)+TOLD(i,j)/4*dy));
end
end
T(:,j)=TriDiag(A,B,C,D); %This 'T' is a function of integral w.r.to y/ how to solve this case usign this 'T'.
dt=0.2+dt;
TOLD=T;
max_difference(j) = max(abs(T(:,j)-T(:,j-1))./max(1,abs(T(:,j-1))));
if max_difference(j)<tol
break
end
max_Iteration=max_difference;
end
end
I want to integrate this 'T' w.r.to Y
T(:,j)=TriDiag(A,B,C,D);
kindly give advice. thank you.
  6 Commenti
Nathan Hardenberg
Nathan Hardenberg il 4 Lug 2023
Sorry I'll not answer the question, since I still do not understand it and your code in the new question was still not executable. Maybe try to reduce your problem. Then other users might be able to help

Accedi per commentare.

Risposte (0)

Community Treasure Hunt

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

Start Hunting!

Translated by