问题遇到的现象和发生背景
问题相关代码,请勿粘贴截图
clear all;clc;close all;
x_min = -1; % 输入最小值,
x_max = 1/4; % 输入最大值,
alpha = @(x,t) t.(2 + cos(2pix));
beta = @(x,t) t.(1 + cos(2pix));
gamma = @(x,t) t.(2 + cos(4pi*x));
Uexact = @(x,t) t.sin(2pix);
Uexact_x = @(x,t) t.2picos(2pix);
Uexact_t = @(x,t) sin(2pix);
Uexact_xx = @(x,t) -t.4pi^2sin(2pi*x);
g1 = @(t) Uexact(x_min,t);
g2 = @(t) Uexact(x_max,t);
fx = @(x,t) Uexact_t(x,t) - alpha(x,t).*Uexact_xx(x,t) - beta(x,t)*Uexact_x(x,t) - gamma(x,t).*Uexact(x,t);
T = 1;
Cons= 0.2;
ArrayN = [40 80 160 320 640]; % 对于dt的收敛;
% AttayN = [10 20 40 80 160]; % 对于dx的收敛使用;
ArrayM = zeros(1,length(ArrayN));
counter = 1;
for N=ArrayN
N = 40;
M = ceil(T/(Cons*(x_max-x_min)^2)*N^2);
% N=80;
% ArrayM(counter)=M;
h = (x_max-x_min)/(N-1); % 步长h
x = zeros(1,N); % 中间点
for i = 1:N
x(i)=x_min+(i-1)*h;
end
dt = T/(M-1); % 步长h
t = zeros(1,M); % 中间点
for n = 1:M
t(n)=(n-1)*dt;
end
u1 = Uexact(x(2:end-1),0).';
A = zeros(N-2,N-2); % 定义矩阵A
A(1,1)=-2;
A(1,2)=1;
A(N-2,N-2)=-2;
A(N-2,N-3)=1;
for j = 2:N-3
A(j,j-1)=1;
A(j,j)=-2;
A(j,j+1)=1;
end
for n = 2:M
tn=(n-1)dt+dt;
% alphau_xx
AlphaMatrix = diag(alpha(x(2:end-1),tn));
Atilda = AlphaMatrix*A;
% beta*u_x
BetaMatrix = diag(beta(x(2:end-1),tn));
C=zeros(N-2,N-2);
C(1,1)=0;
C(1,2)=1;
C(N-2,N-2)=0;
C(N-2,N-3)=-1;
BetaMatrix = diag(beta(x(2:end-1),tn));
for j = 2:N-3
C(j,j+1)=1;
C(j,j)=0;
C(j,j-1)=-1;
end
Ctilda = BetaMatrix*C;
% Gamma*u
GammaMatrix = diag(gamma(x(2:end-1),tn));
D=eye(N-2,N-2);
Dtilda = GammaMatrix*D;
% Implement BCs
b = zeros(N-2,1);
b(1) = (dt/h^2*alpha(x(2),tn) - beta(x(2),tn)*dt/(2*h))*g1(tn);
b(end) = (dt/h^2*alpha(x(end-1),tn) + beta(x(end-1),tn)*dt/(2*h))*g2(tn);
L=D - dt/h^2*Atilda - dt*Dtilda - dt/(2*h)*Ctilda;
u2 = L\(u1 + dt*fx(x(2:end-1),tn) + b);
u1=u2;
Error(counter,n-1) = max(max(abs(Uexact(x(2:end-1),tn).'-u2)));
end
plot(x(2:end-1),u2,'bo-',x(2:end-1),Uexact(x(2:end-1),T),'r*-');
Error2(counter) = max(abs(Error(counter,:)));
counter = counter + 1;
end
loglog(2./ArrayN,Error, 'o--',2./ArrayN,(2./ArrayN).^2,'*--');
xlabel('h')
ylabel('Error')
legend('Exact','Numerical');
grid on;
title('Stability of Implicit Euler method with respect to constant C','interpreter','latex')
loglog((x_max-x_min)./ArrayN,Error2, 'o--',(x_max-x_min)./ArrayN,((x_max-x_min)./ArrayN).^2,'*−−');
xlabel('$dx$','interpreter','latex');
ylabel('Error','interpreter','latex');
title('Convergence with respect to $dx$','interpreter','latex');
legend('$|u(x_i,T)-u_{i,M}|_{\infty}$','$O(dx^2)$','interpreter','latex');
set(gca,'fontsize',13);
grid on;
figure;
loglog(T./ArrayM,Error2,'o--',T./ArrayM,(T./ArrayM),'*--');
xlabel('$dt$','interpreter','latex');
ylabel('Error','interpreter','latex');
title('Convergence with respect to $dt$','interpreter','latex');
legend('$|u(x_i,T)-u_{i,M}|_{\infty}$','$O(dt)$','interpreter','latex');
set(gca,'fontsize',13);
grid on;