upup，我試著寫了一下BTCS的格式，但是不知道怎么加周期性邊界條件[委屈][委屈]，可以幫我看看嘛[呲牙]

隱格式：

時間上采用后向差分，空間上采用中心差分：令CFL=Δt/Δx

格式：u(j,n+1)+CFL/2*[ u(j+1,n+1)- u(j-1,n+1)]= u(j,n)

方程組：A*U(n+1)=U(n)

A=diag(ones(1,Nx+1))+diag(-1/2*CFL*ones(1,Nx),-1)+diag(1/2*CFL*ones(1,Nx),1)

U(n+1)=[u(1,n+1), u(2,n+1)…u(N,n+1)]

U(n)=[u(1,n)+CFL/2*u(0,n+1),u(2,n)…u(N-1,n), u(N,n)- CFL/2*u(N+1,n+1)]

代碼如下

clear

xa = 0; xb = 2*pi; % 求解域

Nx = 200; % 單元數(shù)

Nx1 = Nx + 1; % 格點數(shù)

hx = (xb - xa)/Nx; % Δx

Nt = 500; % 單元數(shù)

tend = 2*pi; % 終止時間

Nt1 = Nt + 1; % 格點數(shù)

ht = (tend - 0)/Nt; % Δx

% 真解

ureal = @(x,t) sin(x-t);?

bcL = 1; bcR = 1;?

CFL = ht/hx; % CFL數(shù)，Δt/Δx

X = xa:hx:xb; % 網(wǎng)格

% 系數(shù)矩陣

U=zeros(Nx1,Nt1);

U(:,1)=ureal(X,0);

A= zeros(Nx1,Nx1);

A = A + diag(ones(1,Nx + 1));

A = A + diag(-1/2*CFL*ones(1,Nx),-1);

A = A + diag(1/2*CFL*ones(1,Nx),1);

F = zeros(Nx1,1); % 右端項

for i=2:Nt1?

??F(1)=U(2,i-1)+1/2*CFL*U(1,i);

??F(end)=U(end-1,i-1)-1/2*CFL*U(end,i);

??for j=2:Nx1-1

????F(j)=U(j,i-1);

??end

??U(:,i)=A\F;

end

figure

plot(X,U(:,end),'ro',X,ureal(X,tend),'b-')

legend('num','exact');