【发布时间】:2014-04-09 20:29:33
【问题描述】:
我不确定是否是这种典型行为,但我正在使用反向差分方法解决有限差分问题。
我用适当的对角项填充了一个稀疏矩阵(沿着中心对角线和它的上方和下方),我尝试使用 MATLAB 的内置方法 (B=A\x) 解决问题,看起来 MATLAB 刚刚得到它错了。
此外,如果使用 inv() 并使用三对角矩阵的逆矩阵,我会得到正确的解决方案。
为什么会这样?
附加信息:
http://pastebin.com/AbuEW6CR(值是标签,所以更容易阅读) 刚度矩阵K:
1 0 0 0
-0.009 1.018 -0.009 0
0 -0.009 1.018 -0.009
0 0 0 1
d 的值:
0
15.55
15.55
86.73
内置输出:
-1.78595556155136e-05
0.00196073713853244
0.00196073713853244
0.0108149483252210
使用 inv(K) 输出:
0
15.42
16.19
86.73
手动输出:
0
15.28
16.18
85.16
代码
nx = 21; %number of spatial steps
nt = 501; %number of time steps (varies between 501 and 4001)
p = alpha * dt / dx^2; %arbitrary constant
a = [0 -p*ones(1,nx-2) 0]'; %diagonal below central diagonal
b = (1+2*p)*ones(nx,1); %central diagonal
c = [1 -p*ones(1,nx-2) 1]'; %diagonal above central diagonal
d = zeros(nx, 1); %rhs values
% Variables a,b,c,d are used for the manual tridiagonal method for
% comparison with MATLAB's built-in functions. The variables represent
% diagonals and the rhs of the matrix
% The equation is K*U(n+1)=U(N)
U = zeros(nx,nt);
% Setting initial conditions
U(:,[1 2]) = (60-32)*5/9;
K = sparse(nx,nx);
% Indices of the sparse matrix which correspond to the diagonal
diagonal = 1:nx+1:nx*nx;
% Populating diagonals
K(diagonal) =1+2*p;
K(diagonal(2:end)-1) =-p;
K(diagonal(1:end-1)+1) =-p;
% Applying dirichlet condition at final spatial step, the temperature is
% derived from a table for predefined values during the calculation
K(end,end-1:end)=[0 1];
% Applying boundary conditions at first spatial step
K(1,1:2) = [1 0];
% Populating rhs values and applying boundary conditions, d=U(n)
d(ivec) = U(ivec,n);
d(nx) = R; %From table
d(1) = 0;
U(:,n+1) = tdm(a,b,c,d); % Manual solver, gives correct answer
U(:,n+1) = d\K; % Built-in solver, gives wrong answer
【问题讨论】:
-
请举一个具体的例子(附代码和结果)。
-
我已经添加了一个,我花了一段时间才把我需要的所有东西放在一起。从我发布的结果来看,使用斜线似乎会产生完全错误的结果。
-
没关系...有史以来最大的脸...我做错了手术...抱歉浪费您的时间。
-
哈哈我正要回答,告诉你你的内置解决方案大致等于 inv() 解决方案除以大约 8000