function shuizhi
E(18)=0;
E(:)=2;
U(18)=0;
U(:)=5;
K=0.0151;
Cj(18)=0;
Cj(1)=10;
Dt=0.1;
Dx=0.5;
j=0;
while (j<10)
Cj=HLSZAD(E,U,Dt,Dx,K,Cj);
j=j+1;
if j>=10
Cj(1)=0;
end
n=length(Cj);
for i=1:n
end
end
function [Cj]=HLSZAD(E,U,Dt,Dx,K,Cj)
%E- j?刻各断面的?向?散系数
%U-j?刻各断面的平均流速
%Dt-????
%Dx-空???
%K-?合衰?系数
%Cj-j?刻各断面?度
n=length(U)-1;
a(n)=0;
b(n)=0;
c(n)=0;
d(n)=0;
for i=1:n
a(i)=-E(i+1)/Dx^2;
b(i)=1/Dt+2*E(i+1)/Dx^2+K/2;
c(i)=-E(i+1)/Dx^2;
d(i)=Cj(i+1)*(1/Dt-U(i+1)/Dx)+Cj(i)*(U(i+1)/Dx-K/2);
disp(d)
end
A=zeros(n,n);
for i=1:n
A(i,i)=b(i);
if(i>1)
A(i,i-1)=a(i);
end
if(i<n)
A(i,i+1)=c(i);
end
end
A(n,n-1)=a(n)-c(n);
A(n,n)=b(n)+2*c(n);
d(1)=d(1)-a(1)*Cj(1);
CK=fzhuigan(A,d);%解三?角方程
Cj(2:n+1)=CK;
function [x]=fzhuigan(A,b)
n=rank(A);
for i=1:n-1
m=A(i+1,i)/A(i,i);
A(i+1,i:i+1)=A(i+1,i:i+1)-m*A(i,i:i+1);
b(i+1)=b(i+1)-m*b(i);
end
x=zeros(1,n);
x(n)=b(n)/A(n,n);
for i=n-1:-1:1
x(i)=(b(i)-A(i,i+1)*x(i+1))/A(i,i);
end
评论0