椭圆型差分

来源：互联网发布：网络直播遗体火化编辑：程序博客网时间：2024/06/10 04:40

椭圆型差分

//精确解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);figuremesh( 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

0 0