椭圆型差分

椭圆型差分       //精确解 function z = gFun(x,y) z = (%e)^( %pi * ( x + y)) * sin( %pi * x) * sin( %pi * y); endfunction //外力 function z = fFun(x,y) z = 2 * %pi * %pi * (%e)^( %pi * ( x + y)) .* ( sin( %pi * x) .* cos( %pi * y) + cos( %pi * x) .* sin( %pi * y)); endfunction //稀疏矩阵 function y = createMatrixA( M) B = diag( ones(M,1) * 4) + diag( -1 * ones(M -1,1), 1) + diag( -1 * ones(M-1,1), -1); diagB = []; for i = 1 : M diagB = sysdiag(diagB,B); end A = diagB; LA = diag( -1 *ones ( M * M - M, 1), -1 * M); UA = diag( -1 * ones( M * M - M ,1), M); A = A + LA + UA ; y = A; endfunction //direclet边界 function y = boundary(x0,y0) y = 0; endfunction //右向量 function y = createVectorB(M, h,X, Y) B = zeros(M * M , 1); for i = 1 : M for j = 1 : M tp = 0; tp = -1 * h * h * fFun(X (i + 1), Y(i + 1)); if ( i == 1) then tp = tp + boundary(X(i + 1 - 1),Y(i + 1)); end if ( i == M) then tp = tp + boundary(X(i + 1 + 1), Y(i + 1)); end if ( j == 1) then tp = tp + boundary(X(i + 1), Y(i + 1 - 1)); end if ( j == M) then tp = tp + boundary(X(i + 1), Y(i + 1 + 1)); end B( (i -1) * M + j) = tp; end end y= B; endfunction //网格剖分 M = 30; X = linspace(0,1, M + 2); Y = linspace(0,1, M + 2); A = createMatrixA(M); //A = sparse(A); b = createVectorB(M, X(2) - X(1), X,Y); solution = A \b; solution = matrix( solution, M, M); globalXY = feval( X(2:M + 1), Y(2:M+1), gFun); mesh( X(2:M+1),Y(2:M+1), globalXY); figure mesh( X(2:M+1), Y(2:M+1), solution); err = abs(solution - globalXY); //figure //mesh(X(2:M+1), Y(2:M+1), err) //plot3d(X(2:M+1), Y(2:M+1),err)   注：1.微分方程差分形式数值解到此为止        2. 椭圆型差分，求解稀疏矩阵，可采用gauss-sidel,jacobi迭代求解，特殊的共轭梯度法；        3, 双曲型守恒性方程、差分收敛性不再累赘。  双曲型 守恒性方程有点深度，有时间再搞；        4. 差分形式对于边界处理繁琐，精度低，但简单易懂；        5. 有限元处理偏微分涉及到图形学相关，funny