reversible matrix operation

4 次查看（过去 30 天）
Peter
Peter 2012-5-16
I am debugging my cfd code and I am checking the validity of my matrices.
I have created a finite difference laplacian matrix in cylindrical coordinates, L. I expect if I take the laplacian of a field, and then multiply by the inverse of the laplacian matrix I should be able to return the exact field that I had before. However, for some reason this is not happening. Basically I am performing the following operations.
L*A = B; C = L\B-A;
I expect that C should equal zero, but it doesn't.
I attached some code showing my problem. Can someone explain what is going wrong?
function [pdiff] = ptest()
AR = 1;
nx = 10;
ny = 20;
lx = 1;
ly = 1;
x = linspace(0,lx,nx+1); hx = lx/nx;
y = linspace(0,ly,ny+1); hy = ly/ny;
[X,Y] = meshgrid(y,x);
pcordr = avg(avg(X')');
pcordz = avg(avg(Y')');
rp = ones(ny,1)*avg(y); rp = rp';
pcenter = pcordr.^2.*pcordz.^2;
plong = reshape(pcenter,[],1);
% % Laplacian matrix
Lp = AR*kron(speye(ny),K1(nx,hx,1))+(1/AR)*kron(K1(ny,hy,1),speye(nx))+...
(1/AR)*kron((1./rp).*Kp1(ny,hy),speye(nx));
% Test for reverse operations
pfield = Lp*plong;
pdiff = Lp\pfield - plong;
function A = K1(n,h,a11)
% a11: Neumann=1, Dirichlet=2, Dirichlet mid=3;
A = spdiags([-1 a11 0;ones(n-2,1)*[-1 2 -1];0 a11 -1],-1:1,n,n)'/h^2;
function A = Kp1(n,h)
% a11: Neumann=1, Dirichlet=2, Dirichlet mid=3;
A = spdiags([-1 1 0;ones(n-2,1)*[-1 0 1];0 -1 1],-1:1,n,n)'/(h*2);
function B = avg(A,k)
if nargin<2, k = 1; end
if size(A,1)==1, A = A'; end
if k<2, B = (A(2:end,:)+A(1:end-1,:))/2; else, B = avg(A,k-1); end
if size(A,2)==1, B = B'; end

请先登录，再回答此问题。

回答（3 个）

Walter Roberson
Walter Roberson 2012-5-16
The output that results, C: is it on the order of 1e-16 times the largest (absolute value) input? If so then you are observing round-off error.
  2 个评论
Peter
Peter 2012-5-16
no, its considerably larger
Walter Roberson
Walter Roberson 2012-5-16
In my experimentation, the largest I have found has been less than the largest of max(eps(A*x)), max(eps(A)), max(eps(x))

请先登录，再进行评论。

Richard Brown
Richard Brown 2012-5-16
Your code doesn't run for me -- first, I don't have a function avg - did you mean to use mean?
Second, if I replace all occurences of avg with mean, the line defining the Laplacian doesn't work - in particular you can't do
(1./rp).*Kp1(ny,hy)
as your function Kp1 defines a matrix and 1./rp is a vector ...
Richard Brown
Richard Brown 2012-5-16
Your Laplacian matrix is not full rank. Because your system is consistent, you have multiple solutions, and you're just getting one of them when you do Lp \ pfield
Unfortunately because your matrix Lp is sparse, you don't see the usual warning message about the matrix being badly scaled or singular ...

类别

在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by