function [t,f0,f]=flplot(dp,opt)
n=20;
h=1/n;
for i=1:n+1
t(i)=h*(i-1);
f0(i)=(exp(dp*(i-1))/(exp(dp*n)));
end
q0=100;
qL=1000;
switch num2str(opt)
case'1'
%中心差分
aE=1-0.5*dp;
aW=1+0.5*dp;
case'2'
%一阶迎风格式
aE=1;
aW=1+dp;
case'3'
%混合格式
aE=max(0,(1-0.5*dp));
aW=max(0,(1-0.5*dp))+dp;
case'4'
%quick
aE=1-(3/8)*dp;
aW=1+(7/8)*dp;
aWW=-(1/8)*dp;
end
A=zeros(n-1);
if opt<4
aP=aE+aW;
for i=1:n-1
A(i,i)=-aP;
if(i>1)
A(i,i-1)=aW;
end
if(i<n-1)
A(i,i+1)=aE;
end
end
b=zeros(n-1,1);
b(1,1)=-q0*aW;
b(n-1,1)=-qL*aE;
else
aP=aWW+aW+aE;
A(1,1)=-(2+dp);
A(1,2)=1;%一阶迎风格式
for i=2:n-1
A(i,i)=-aP;
A(i,i-1)=aW;
if(i>2)
A(i,i-2)=aWW;
end
if(i<n-1)
A(i,i+1)=aE;
end
end
b=zeros(n-1,1);
b(1,1)=-q0*(1+dp);
b(2,1)=-q0*aWW;
b(n-1,1)=-qL*aE;
end
q=A\b;
f=zeros(n+1,1);
for i=1:n-1
f(i+1)=(q(i)-q0)/(qL-q0);
end
f(n+1)=1;
评论2