有一个n×m的网格,小A在(x, y), 有另外k个人,第i个人在(xi, yi), 每一秒每个人必须移动到相邻的格子,判断是否存在一个时刻小A和至少一个人在同一个格子。

这段C++代码是用于解决一个抓人游戏问题,判断第i个人能否抓住小A,条件是其与小A在x轴和y轴上的绝对差之和为偶数。如果满足这个条件的人数超过1,表示无法抓住小A;否则可以抓住。

题目

思路:独立地看第i个人,当且仅当(abs(x[i] - x) + abs(y[i] - y)) % 2 == 0, 第i个人才能抓到小A

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define pb push_back
#define fi first
#define se second
#define lson p << 1
#define rson p << 1 | 1
const int maxn = 1e6 + 5, inf = 1e18, maxm = 4e4 + 5, mod = 1e9 + 7, N = 1e6;
int a[maxn], b[maxn];
// bool vis[maxn];
int n, m;
int vis[105][105] = {0};
string s;

void solve(){
    int res = 0;
    int q, k;
    int x0, y0;
    int x[105], y[105];
    cin >> n >> m >> k;
    cin >> x[0] >> y[0];
    int cnt = 0;
    for(int i = 1; i <= k; i++){
        cin >> x[i] >> y[i];
        int dx = abs(x[0] - x[i]);
        int dy = abs(y[0] - y[i]);
        if((dx + dy) % 2 == 0){
            cnt++;
        }
    }
    if(cnt >= 1){
        cout << "No\n";
    }
    else cout << "Yes\n";
    
}
    
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0);
    
    int T = 1;
    cin >> T;
    while (T--)
    {
        solve();
    }
    return 0;
}

``` #include<bits/stdc++.h> using namespace std; #define ll long long const int N=2000+5; struct lianb{ int n; int a[N]; int ne[N]; int pr[N]; lianb(){ for(int i=1;i<N;++i){ a[i]=ne[i]=pr[i]=0; } } int pre(int x,int y){ y-=a[x]; while(y>0){ x=pr[x]; y-=a[x]; } return x; } int nex(int x,int y){ y-=a[x]; while(y>0){ x=ne[x]; y-=a[x]; } return x; } int pr_len(int x){ return x-pr[x]; } int ne_len(int x){ return ne[x]-x; } int run(int wn,int k){ n=wn; int la=0; pr[0]=0; a[0]=a[n+1]=1e9; for(int i=1;i<=n;++i){ if(a[i]!=0){ ne[la]=i; pr[i]=la; la=i; } } ne[la]=ne[n+1]=n+1; int ans=0; for(int l=ne[0];l<=n;l=ne[l]){ ans+=pr_len(l)*(n-nex(l,k)+1); } return ans; } }; int n,m,w,k; vector<int>v[N]; lianb ls; int main(){ ios::sync_with_stdio(0); cin.tie(0); cout.tie(0); cin>>n>>m>>w>>k; for(int i=1;i<=w;++i){ int x,y; cin>>x>>y; v[x].push_back(y); } ll ans=0; for(int i=1;i<=n;++i){ ls=lianb(); for(int j=i;j<=n;++j){ for(int x:v[j]){ ++ls.a[x]; } } int wp=ls.run(m,k); for(int j=n;j>=i;--j){ ans+=wp; for(int x:v[j]){ for(int l=ls.pre(x,k);l<=x;l=ls.ne[l]){ wp-=ls.pr_len(l)*ls.ne_len(ls.nex(l,k)); } } for(int x:v[j]){ --ls.a[x]; } } } cout<<ans<<endl; return 0; }```debug 保罗在交响乐团。弦乐部分排列成一个 r × c 矩形网格,除了 n 位中提琴手外,其他都是小提琴手。保罗非常喜欢中提琴,所以他想拍一张照片,至少要包括 k 位中提琴手。保罗可以拍摄交响乐团中的任何轴平行矩形。请计算保罗可以拍摄的可能照片数量。 如果对应矩形的坐标不同,则认为两张照片是不同的。 输入 输入的第一行包含四个用空格分隔的整数 r, c, n, k (1 ≤ r, c, n ≤ 3000, 1 ≤ k ≤ min(n, 10)) — 弦乐部分的行数列数,中提琴的总数,保罗想要照片中至少包括的中提琴的数量。 接下来的 n 行,每行包含两个整数 xi yi (1 ≤ xi ≤ r, 1 ≤ yi ≤ c):第 i 位中提琴的位置。保证输入中没有重复的位置。 输出 输出一个整数 — 保罗可以拍摄的照片数量,其中至少包括 k 位中提琴手。 Verdict: WRONG_ANSWER Input 10 10 10 4 7 9 8 1 5 7 6 8 6 4 5 8 4 3 1 3 5 1 6 9 Participant's output 598 Jury's answer 619 Checker comment wrong answer 1st numbers differ - expected: '619', found: '598'
04-04
tic; clear clc % 追赶法函数定义 function x = chase(a, b, c, d) n = length(b); % 前向消元 for i = 2:n m = a(i-1) / b(i-1); b(i) = b(i) - m * c(i-1); d(i) = d(i) - m * d(i-1); end % 回代求解 x = zeros(n, 1); x(n) = d(n) / b(n); for i = n-1:-1:1 x(i) = (d(i) - c(i) * x(i+1)) / b(i); end end M=[10,20,40,80];% 空间网格数量 N=M;% 时间网格数量 for p=1:length(M) h=1/M(p);% 这里定义空间步长等距 tau=1/N(p); % 时间步长 x=0:h:1; y=0:h:1; t=0:tau:1; Numerical(M(p)+1,M(p)+1,N(p)+1)=0;%u numerical(M(p)+1,M(p)-1)=0;%u* % 求解u*ijuij过程中构造三对角矩阵 % a 表示下对角线元素 % b 表示主对角线元素 % c 表示上对角线元素 a=-tau^2/(2*h^2)*ones(M(p)-2,1); b=(tau^2/h^2+1)*ones(M(p)-1,1); c=-tau^2/(2*h^2)*ones(M(p)-2,1); for i=1:M(p)+1 for j=1:M(p)+1 Numerical(i,j,1)=exp(1/2*(x(i)+y(j)));% 初值Numerical(x,y,0)=u(i,j,0) end end for j=1:M(p)+1 for k=1:N(p)+1 Numerical(1,j,k)=exp(1/2*y(j)-t(k));% 边值Numerical(0,y,t)=u(0,j,k) end end for j=1:M(p)+1 for k=1:N(p)+1 Numerical(M(p)+1,j,k)=exp(1/2*(1+y(j))-t(k));% 边值Numerical(1,y,t)=u(m,j,k) end end for i=1:M(p)+1 for k=1:N(p)+1 Numerical(i,1,k)=exp(1/2*x(i)-t(k));% 边值Numerical(x,0,t)=u(i,0,k) end end for i=1:M(p)+1 for k=1:N(p)+1 Numerical(i,M(p)+1,k)=exp(1/2*(1+x(i))-t(k));% 边值Numerical(x,1,t)=u(i,m,k) end end % f=inline('-1/2*exp(1/2*(x+y)-t)','x','y','t');% f(x,y,t) % fun=inline('exp(1/2*(x+y)-t)','x','y','t');% 精确解 f = @(x, y, t) 0.5 * exp(0.5*(x+y) - t); % 源项 fun = @(x, y, t) exp(0.5*(x+y) - t); % 精确解 psi = @(x,y) -exp(1/2*(x+y)); rou = @(x,y) -exp(1/2*(x+y)); % 计算精确解 for i=1:M(p)+1 for j=1:M(p)+1 for k=1:N(p)+1 Accurate(i,j,k)=fun(x(i),y(j),t(k)); end end end % 核心部分(第一层) for k=1:N(p) for j=1:M(p)-1% 固定j % 动态计算 u* 的边界 numerical(1,j)=-tau^2/(2*h^2)*Numerical(1,j,2)+(tau^2/h^2+1)*Numerical(1,j+1,2)-tau^2/(2*h^2)*Numerical(1,j+2,2);% u*0j numerical(M(p)+1,j)=-tau^2/(2*h^2)*Numerical(M(p)+1,j,2)+(tau^2/h^2+1)*Numerical(M(p)+1,j+1,2)-tau^2/(2*h^2)*Numerical(M(p)+1,j+2,2);% u*mj % 循环生成1式右端列向量 for i=1:M(p)-1 numerical_right_vector(i,1)=Numerical(i+1,j+1,1)+tau*psi(x(i+1),y(j+1))-tau^3*rou(x(i+1),y(j+1))/3+... tau^2*f(x(i+1),y(j+1),t(2))/2; end % 添加边界贡献 numerical_right_vector(1,1)=numerical_right_vector(1,1)+tau^2/(2*h^2)*numerical(1,j); numerical_right_vector(M(p)-1,1)=numerical_right_vector(M(p)-1,1)+tau^2/(2*h^2)*numerical(M(p)+1,j); numerical(2:M(p),j)=chase(a,b,c,numerical_right_vector); end for i=1:M(p)-1 % 固定i for j=1:M(p)-1 Numerical_right_vector(j,1)=numerical(i+1,j); end Numerical_right_vector(1,1)=Numerical_right_vector(1,1)+tau^2/(2*h^2)*Numerical(i+1,1,2); Numerical_right_vector(M(p)-1,1)=Numerical_right_vector(M(p)-1,1)+tau^2/(2*h^2)*Numerical(i+1,M(p)+1,2); Numerical(i+1,2:M(p),k+1)=chase(a,b,c,Numerical_right_vector); end end %第二层 for k=2:N(p) for j=1:M(p)-1% 固定j numerical(1,j)=-tau^2/(2*h^2)*(Numerical(1,j,k+1)+Numerical(1,j,k-1))/2+... (tau^2/h^2+1)*(Numerical(1,j+1,k+1)+Numerical(1,j+1,k-1))/2+... -tau^2/(2*h^2)*(Numerical(1,j+2,k+1)+Numerical(1,j+2,k-1))/2;% u*0j numerical(M(p)+1,j)=-tau^2/(2*h^2)*(Numerical(M(p)+1,j,k+1)+Numerical(M(p)+1,j,k-1))/2+... (tau^2/h^2+1)*(Numerical(M(p)+1,j+1,k+1)+Numerical(M(p)+1,j+1,k-1))/2+... -tau^2/(2*h^2)*(Numerical(M(p)+1,j+2,k+1)+Numerical(M(p)+1,j+2,k-1))/2;% u*mj % 循环生成1式右端列向量 numerical_right_vector = zeros(M(p)-1, 1); for i=1:M(p)-1 numerical_right_vector(i,1)=Numerical(i+1,j+1,k)+tau^2/2*f(x(i+1),y(j+1),t(k)); end % 添加边界贡献 numerical_right_vector(1,1)=numerical_right_vector(1,1)+tau^2/(2*h^2)*numerical(1,j); numerical_right_vector(M(p)-1,1)=numerical_right_vector(M(p)-1,1)+tau^2/(2*h^2)*numerical(M(p)+1,j); numerical(2:M(p),j)=chase(a,b,c,numerical_right_vector); end for i=1:M(p)-1 % 固定i boundary_left = (Accurate(i+1,1,k+1) + Accurate(i+1,1,k-1))/2; boundary_right = (Accurate(i+1,M(p)+1,k+1) + Accurate(i+1,M(p)+1,k-1))/2; Numerical_right_vector = zeros(M(p)-1, 1); for j=1:M(p)-1 Numerical_right_vector(j,1)=numerical(i+1,j); end Numerical_right_vector(1,1)=Numerical_right_vector(1,1)+tau^2/(2*h^2)*boundary_left; Numerical_right_vector(M(p)-1,1)=Numerical_right_vector(M(p)-1,1)+tau^2/(2*h^2)*boundary_right; temp_result=chase(a,b,c,Numerical_right_vector); Numerical(i+1,2:M(p),k+1) = 2*temp_result' - Numerical(i+1,2:M(p),k-1); end end error=abs(Numerical(:,:,M(p)+1)-Accurate(:,:,M(p)+1)); error_inf(p)=max(max(error)); fprintf('网格 %dx% d 计算完成,最大误差: %.6e\n', M(p), M(p), error_inf(p)); figure(p) [X,Y]=meshgrid(y,x); subplot(1,3,1) surf(X,Y,Accurate(:,:,end)); xlabel('x');ylabel('y');zlabel('u'); title('精确解图像'); grid on; subplot(1,3,2) surf(X,Y,Numerical(:,:,end)); xlabel('x');ylabel('y');zlabel('u'); title('数值解图像'); grid on; subplot(1,3,3) surf(X,Y,error); xlabel('x');ylabel('y');zlabel('error'); title('误差图像'); grid on; sgtitle(sprintf('D''Yakonov ADI 方法结果 (网格: %dx%d)', M(p), M(p))); end % 计算收敛阶 fprintf('\n=== 收敛性分析 ===\n'); fprintf('%-10s %-12s %-15s\n', '网格尺寸', '最大误差', '收敛阶'); fprintf('%s\n', repmat('-', 1, 40)); fprintf('%-10d %-12.2e %-15s\n', M(1), error_inf(1), '--'); for k=2:length(M) X=error_inf(k-1)/error_inf(k); Norm(k-1)=X; fprintf('%-10d %-12.2e %-15.4f\n', M(k), error_inf(k), Norm(k-1)); end figure(length(M)+1) plot(1:length(Norm),Norm,'-b^','LineWidth',2,'MarkerSize',8); xlabel('网格加密序列');ylabel('误差阶数'); title('D''Yakonov ADI格式误差阶'); grid on; % 在图上标注数值 for i=1:length(Norm) text(i, Norm(i), sprintf('%.3f', Norm(i)), ... 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'center'); end fprintf('\n=== 计算总结 ===\n'); fprintf('✅ D''Yakonov ADI 方法成功实现(使用追赶法)\n'); fprintf('✅ 在计算过程中动态计算 u* 的边界\n'); fprintf('✅ 所有时间层空间点的数值解、精确解误差已计算\n'); fprintf('✅ 最大误差分析完成,收敛性验证通过\n'); toc;将其改为julia语言
最新发布
10-27
tic; clear clc % 追赶法函数定义 function x = chase(a, b, c, d) n = length(b); % 前向消元 for i = 2:n m = a(i-1) / b(i-1); b(i) = b(i) - m * c(i-1); d(i) = d(i) - m * d(i-1); end % 回代求解 x = zeros(n, 1); x(n) = d(n) / b(n); for i = n-1:-1:1 x(i) = (d(i) - c(i) * x(i+1)) / b(i); end end M=[10,20,40,80];% 空间网格数量 N=M;% 时间网格数量 for p=1:length(M) h=1/M(p);% 这里定义空间步长等距 tau=1/N(p); % 时间步长 x=0:h:1; y=0:h:1; t=0:tau:1; Numerical(M(p)+1,M(p)+1,N(p)+1)=0;%u numerical(M(p)+1,M(p)-1)=0;%u* % 求解u*ijuij过程中构造三对角矩阵 % a 表示下对角线元素 % b 表示主对角线元素 % c 表示上对角线元素 a=-tau/(2*h^2)*ones(M(p)-2,1); b=(tau/h^2+1)*ones(M(p)-1,1); c=-tau/(2*h^2)*ones(M(p)-2,1); for i=1:M(p)+1 for j=1:M(p)+1 Numerical(i,j,1)=exp(1/2*(x(i)+y(j)));% 初值Numerical(x,y,0)=u(i,j,0) end end for j=1:M(p)+1 for k=1:N(p)+1 Numerical(1,j,k)=exp(1/2*y(j)-t(k));% 边值Numerical(0,y,t)=u(0,j,k) end end for j=1:M(p)+1 for k=1:N(p)+1 Numerical(M(p)+1,j,k)=exp(1/2*(1+y(j))-t(k));% 边值Numerical(1,y,t)=u(m,j,k) end end for i=1:M(p)+1 for k=1:N(p)+1 Numerical(i,1,k)=exp(1/2*x(i)-t(k));% 边值Numerical(x,0,t)=u(i,0,k) end end for i=1:M(p)+1 for k=1:N(p)+1 Numerical(i,M(p)+1,k)=exp(1/2*(1+x(i))-t(k));% 边值Numerical(x,1,t)=u(i,m,k) end end % f=inline('-3/2*exp(1/2*(x+y)-t)','x','y','t');% f(x,y,t) % fun=inline('exp(1/2*(x+y)-t)','x','y','t');% 精确解 f = @(x, y, t) -1.5 * exp(0.5*(x+y) - t); % 源项 fun = @(x, y, t) exp(0.5*(x+y) - t); % 精确解 % 计算精确解 for i=1:M(p)+1 for j=1:M(p)+1 for k=1:N(p)+1 Accurate(i,j,k)=fun(x(i),y(j),t(k)); end end end % 核心部分 for k=1:N(p) for j=1:M(p)-1% 固定j % 动态计算 u* 的边界 numerical(1,j)=-tau/(2*h^2)*Numerical(1,j,k+1)+(tau/h^2+1)*Numerical(1,j+1,k+1)-tau/(2*h^2)*Numerical(1,j+2,k+1);% u*0j numerical(M(p)+1,j)=-tau/(2*h^2)*Numerical(M(p)+1,j,k+1)+(tau/h^2+1)*Numerical(M(p)+1,j+1,k+1)-tau/(2*h^2)*Numerical(M(p)+1,j+2,k+1);% u*mj % 循环生成1式右端列向量 for i=1:M(p)-1 numerical_right_vector(i,1)=tau*f(x(i+1),y(j+1),t(k)+tau/2)+Numerical(i+1,j+1,k)... +tau/(2*h^2)*(Numerical(i,j+1,k)-2*Numerical(i+1,j+1,k)+Numerical(i+2,j+1,k))... +tau/(2*h^2)*(Numerical(i+1,j,k)-2*Numerical(i+1,j+1,k)+Numerical(i+1,j+2,k))... +tau^2/(4*h^4)*(Numerical(i,j,k)-2*Numerical(i+1,j,k)+Numerical(i+2,j,k))... +tau^2/(4*h^4)*(-2*Numerical(i,j+1,k)+4*Numerical(i+1,j+1,k)-2*Numerical(i+2,j+1,k))... +tau^2/(4*h^4)*(Numerical(i,j+2,k)-2*Numerical(i+1,j+2,k)+Numerical(i+2,j+2,k)); end % 添加边界贡献 numerical_right_vector(1,1)=numerical_right_vector(1,1)+tau/(2*h^2)*numerical(1,j); numerical_right_vector(M(p)-1,1)=numerical_right_vector(M(p)-1,1)+tau/(2*h^2)*numerical(M(p)+1,j); numerical(2:M(p),j)=chase(a,b,c,numerical_right_vector); end for i=1:M(p)-1 % 固定i for j=1:M(p)-1 Numerical_right_vector(j,1)=numerical(i+1,j); end Numerical_right_vector(1,1)=Numerical_right_vector(1,1)+tau/(2*h^2)*Numerical(i+1,1,k+1); Numerical_right_vector(M(p)-1,1)=Numerical_right_vector(M(p)-1,1)+tau/(2*h^2)*Numerical(i+1,M(p)+1,k+1); Numerical(i+1,2:M(p),k+1)=chase(a,b,c,Numerical_right_vector); end end error=abs(Numerical(:,:,end)-Accurate(:,:,end)); error_inf(p)=max(max(error)); fprintf('网格 %dx%d 计算完成,最大误差: %.6e\n', M(p), M(p), error_inf(p)); figure(p) [X,Y]=meshgrid(y,x); subplot(1,3,1) surf(X,Y,Accurate(:,:,end)); xlabel('x');ylabel('y');zlabel('u'); title('精确解图像'); grid on; subplot(1,3,2) surf(X,Y,Numerical(:,:,end)); xlabel('x');ylabel('y');zlabel('u'); title('数值解图像'); grid on; subplot(1,3,3) surf(X,Y,error); xlabel('x');ylabel('y');zlabel('error'); title('误差图像'); grid on; sgtitle(sprintf('D''Yakonov ADI 方法结果 (网格: %dx%d)', M(p), M(p))); end % 计算收敛阶 fprintf('\n=== 收敛性分析 ===\n'); fprintf('%-10s %-12s %-15s\n', '网格尺寸', '最大误差', '收敛阶'); fprintf('%s\n', repmat('-', 1, 40)); fprintf('%-10d %-12.2e %-15s\n', M(1), error_inf(1), '--'); for k=2:length(M) X=error_inf(k-1)/error_inf(k); Norm(k-1)=X; fprintf('%-10d %-12.2e %-15.4f\n', M(k), error_inf(k), Norm(k-1)); end figure(length(M)+1) plot(1:length(Norm),Norm,'-b^','LineWidth',2,'MarkerSize',8); xlabel('网格加密序列');ylabel('误差阶数'); title('D''Yakonov ADI格式误差阶'); grid on; % 在图上标注数值 for i=1:length(Norm) text(i, Norm(i), sprintf('%.3f', Norm(i)), ... 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'center'); end fprintf('\n=== 计算总结 ===\n'); fprintf('✅ D''Yakonov ADI 方法成功实现(使用追赶法)\n'); fprintf('✅ 在计算过程中动态计算 u* 的边界\n'); fprintf('✅ 所有时间层空间点的数值解、精确解误差已计算\n'); fprintf('✅ 最大误差分析完成,收敛性验证通过\n'); toc改为紧致交替方向隐格式
10-14
tic; clear clc % 追赶法函数定义 function x = chase(a, b, c, d) n = length(b); % 前向消元 for i = 2:n m = a(i-1) / b(i-1); b(i) = b(i) - m * c(i-1); d(i) = d(i) - m * d(i-1); end % 回代求解 x = zeros(n, 1); x(n) = d(n) / b(n); for i = n-1:-1:1 x(i) = (d(i) - c(i) * x(i+1)) / b(i); end end M=[10,20,40,80];% 空间网格数量 N=M;% 时间网格数量 for p=1:length(M) h=1/M(p);% 这里定义空间步长等距 tau=1/N(p); % 时间步长 x=0:h:1; y=0:h:1; t=0:tau:1; Numerical(M(p)+1,M(p)+1,N(p)+1)=0;%u numerical(M(p)+1,M(p)-1)=0;%u* % 求解u*ijuij过程中构造三对角矩阵 % a 表示下对角线元素 % b 表示主对角线元素 % c 表示上对角线元素 a=-tau^2/(2*h^2)*ones(M(p)-2,1); b=(tau^2/h^2+1)*ones(M(p)-1,1); c=-tau^2/(2*h^2)*ones(M(p)-2,1); for i=1:M(p)+1 for j=1:M(p)+1 Numerical(i,j,1)=exp(1/2*(x(i)+y(j)));% 初值Numerical(x,y,0)=u(i,j,0) end end for j=1:M(p)+1 for k=1:N(p)+1 Numerical(1,j,k)=exp(1/2*y(j)-t(k));% 边值Numerical(0,y,t)=u(0,j,k) end end for j=1:M(p)+1 for k=1:N(p)+1 Numerical(M(p)+1,j,k)=exp(1/2*(1+y(j))-t(k));% 边值Numerical(1,y,t)=u(m,j,k) end end for i=1:M(p)+1 for k=1:N(p)+1 Numerical(i,1,k)=exp(1/2*x(i)-t(k));% 边值Numerical(x,0,t)=u(i,0,k) end end for i=1:M(p)+1 for k=1:N(p)+1 Numerical(i,M(p)+1,k)=exp(1/2*(1+x(i))-t(k));% 边值Numerical(x,1,t)=u(i,m,k) end end % f=inline('-1/2*exp(1/2*(x+y)-t)','x','y','t');% f(x,y,t) % fun=inline('exp(1/2*(x+y)-t)','x','y','t');% 精确解 f = @(x, y, t) -0.5 * exp(0.5*(x+y) - t); % 源项 fun = @(x, y, t) exp(0.5*(x+y) - t); % 精确解 psi = @(x,y) -exp(1/2*(x+y)); rou = @(x,y) -exp(1/2*(x+y)); % 计算精确解 for i=1:M(p)+1 for j=1:M(p)+1 for k=1:N(p)+1 Accurate(i,j,k)=fun(x(i),y(j),t(k)); end end end % 核心部分 for k=1:N(p) for j=1:M(p)-1% 固定j % 动态计算 u* 的边界 numerical(1,j)=-tau^2/(2*h^2)*Numerical(1,j,2)+(tau^2/h^2+1)*Numerical(1,j+1,2)-tau^2/(2*h^2)*Numerical(1,j+2,2);% u*0j numerical(M(p)+1,j)=-tau^2/(2*h^2)*Numerical(M(p)+1,j,2)+(tau^2/h^2+1)*Numerical(M(p)+1,j+1,2)-tau^2/(2*h^2)*Numerical(M(p)+1,j+2,2);% u*mj % 循环生成1式右端列向量 for i=1:M(p)-1 numerical_right_vector(i,1)=Numerical(i+1,j+1,1)+tau*psi(x(i+1),y(j+1))-tau^3*rou(x(i+1),y(j+1))/3+... tau^2*f(x(i+1),y(j+1),t(2))/2; end % 添加边界贡献 numerical_right_vector(1,1)=numerical_right_vector(1,1)+tau^2/(2*h^2)*numerical(1,j); numerical_right_vector(M(p)-1,1)=numerical_right_vector(M(p)-1,1)+tau^2/(2*h^2)*numerical(M(p)+1,j); numerical(2:M(p),j)=chase(a,b,c,numerical_right_vector); end for i=1:M(p)-1 % 固定i for j=1:M(p)-1 Numerical_right_vector(j,1)=numerical(i+1,j); end Numerical_right_vector(1,1)=Numerical_right_vector(1,1)+tau^2/(2*h^2)*Numerical(i+1,1,2); Numerical_right_vector(M(p)-1,1)=Numerical_right_vector(M(p)-1,1)+tau^2/(2*h^2)*Numerical(i+1,M(p)+1,2); Numerical(i+1,2:M(p),k+1)=chase(a,b,c,Numerical_right_vector); end end for k=2:N(p) for j=1:M(p)-1% 固定j numerical(1,j)=-tau^2/(2*h^2)*Numerical(1,j,k+1)+... (tau^2/h^2+1)*Numerical(1,j+1,k+1)+... -tau^2/(2*h^2)*Numerical(1,j+2,k+1);% u*0j numerical(M(p)+1,j)=-tau^2/(2*h^2)*Numerical(M(p)+1,j,k+1)+... (tau^2/h^2+1)*Numerical(M(p)+1,j+1,k+1)+... -tau^2/(2*h^2)*Numerical(M(p)+1,j+2,k+1);% u*mj % 循环生成1式右端列向量 for i=1:M(p)-1 numerical_right_vector(i,1)=2*Numerical(i+1,j+1,k)+f(x(i+1),y(j+1),t(k))*tau^2+... tau^2/(2*h^2)*(Numerical(i,j+1,k-1)+Numerical(i+2,j+1,k-1))+... tau^2/(2*h^2)*(Numerical(i+1,j,k-1)+Numerical(i+1,j+2,k-1))-(1+tau^2/(h^2)+tau^2/(h^2))*Numerical(i+1,j+1,k-1); end % 添加边界贡献 numerical_right_vector(1,1)=numerical_right_vector(1,1)+tau^2/(2*h^2)*numerical(1,j); numerical_right_vector(M(p)-1,1)=numerical_right_vector(M(p)-1,1)+tau^2/(2*h^2)*numerical(M(p)+1,j); numerical(2:M(p),j)=chase(a,b,c,numerical_right_vector); end for i=1:M(p)-1 % 固定i for j=1:M(p)-1 Numerical_right_vector(j,1)=numerical(i+1,j); end Numerical_right_vector(1,1)=Numerical_right_vector(1,1)+tau^2/(2*h^2)*Numerical(i+1,1,k+1); Numerical_right_vector(M(p)-1,1)=Numerical_right_vector(M(p)-1,1)+tau^2/(2*h^2)*Numerical(i+1,M(p)+1,k+1); Numerical(i+1,2:M(p),k+1)=chase(a,b,c,Numerical_right_vector); end end error=abs(Numerical(:,:,end)-Accurate(:,:,end)); error_inf(p)=max(max(error)); fprintf('网格 %dx%d 计算完成,最大误差: %.6e\n', M(p), M(p), error_inf(p)); figure(p) [X,Y]=meshgrid(y,x); subplot(1,3,1) surf(X,Y,Accurate(:,:,end)); xlabel('x');ylabel('y');zlabel('u'); title('精确解图像'); grid on; subplot(1,3,2) surf(X,Y,Numerical(:,:,end)); xlabel('x');ylabel('y');zlabel('u'); title('数值解图像'); grid on; subplot(1,3,3) surf(X,Y,error); xlabel('x');ylabel('y');zlabel('error'); title('误差图像'); grid on; sgtitle(sprintf('D''Yakonov ADI 方法结果 (网格: %dx%d)', M(p), M(p))); end % 计算收敛阶 fprintf('\n=== 收敛性分析 ===\n'); fprintf('%-10s %-12s %-15s\n', '网格尺寸', '最大误差', '收敛阶'); fprintf('%s\n', repmat('-', 1, 40)); fprintf('%-10d %-12.2e %-15s\n', M(1), error_inf(1), '--'); for k=2:length(M) X=error_inf(k-1)/error_inf(k); Norm(k-1)=X; fprintf('%-10d %-12.2e %-15.4f\n', M(k), error_inf(k), Norm(k-1)); end figure(length(M)+1) plot(1:length(Norm),Norm,'-b^','LineWidth',2,'MarkerSize',8); xlabel('网格加密序列');ylabel('误差阶数'); title('D''Yakonov ADI格式误差阶'); grid on; % 在图上标注数值 for i=1:length(Norm) text(i, Norm(i), sprintf('%.3f', Norm(i)), ... 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'center'); end fprintf('\n=== 计算总结 ===\n'); fprintf('✅ D''Yakonov ADI 方法成功实现(使用追赶法)\n'); fprintf('✅ 在计算过程中动态计算 u* 的边界\n'); fprintf('✅ 所有时间层空间点的数值解、精确解误差已计算\n'); fprintf('✅ 最大误差分析完成,收敛性验证通过\n'); toc;
10-23
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

__night_

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值