moon 發問時間： 電腦與網際網路軟體 · 1 0 年前

# 請問一維向量、二維矩陣偏微分matlab之語法?

### 1 個解答

• 老師
Lv 7
1 0 年前
最佳解答

給你一個例題參考function [u,x,y,error]=poissonfd(a,c,b,d,nx,ny,f,g,uex)%POISSONFD two-dimensional Poisson solver%   [U,X,Y]=POISSONFD(A,C,B,D,NX,NY,FUN,BOUND) solves with%   the 5-points finite difference scheme the problem%   -LAPL(U) = FUN in the rectangle (A,C)X(B,D) with%   Dirichlet boundary conditions U(X,Y)=BOUND(X,Y) for any%   (X,Y) on the boundary on the rectangle. %%   [U,X,Y,ERROR]=POISSONFD(A,C,B,D,NX,NY,FUN,BOUND,UEX) computes%   also the maximum nodal error ERROR with respect to the exact%   solution UEX.%   FUN, BOUND and UEX can be inline functions.% refer to Alfio Quarteroni, Fausto Saleri,李敏波譯% MATLAB 科學計算(Scientific Computing with MATLAB) 清華大學出版社% p. 174, 程序 18-poissonfd: % 利用五點有限差分法求解含有Dirichlet條件的Poisson問題的近似解if nargin == 8 & nargout==4    uex = '0 + 0.*x'; end nx=nx+1;         ny=ny+1;hx = (b-a)/nx;   hy = (d-c)/ny; nx1 = nx+1;      hx2 = hx^2; hy2 = hy^2; kii = 2/hx2+2/hy2;    kix = -1/hx2;    kiy = -1/hy2;dim = (nx+1)*(ny+1); K = speye(dim,dim); rhs = zeros(dim,1); y = c;for m = 2:ny    x = a; y = y + hy;    for n = 2:nx        i = n+(m-1)*(nx+1);         x = x + hx;         rhs(i) = feval(f,x,y);          K(i,i) = kii;         K(i,i-1) = kix;         K(i,i+1) = kix;        K(i,i+nx1) = kiy;         K(i,i-nx1) = kiy;    endend rhs1 = zeros(dim,1);x = [a:hx:b]; rhs1(1:nx1) = feval(g,x,c);rhs1(dim-nx:dim) = feval(g,x,d);y = [c:hy:d]; rhs1(1:nx1:dim-nx) = feval(g,a,y); rhs1(nx1:nx1:dim) = feval(g,b,y);rhs = rhs - K*rhs1; nbound = [[1:nx1],[dim-nx:dim],[1:nx1:dim-nx],[nx1:nx1:dim]];ninternal = setdiff([1:dim],nbound);K = K(ninternal,ninternal); rhs = rhs(ninternal); utemp = K\rhs; uh = rhs1; uh (ninternal) = utemp; k = 1; y = c;for j = 1:ny+1    x = a;    for i = 1:nx1        u(i,j) = uh(k);         k = k + 1;         ue(i,j) = feval(uex,x,y);         x = x + hx;     end    y = y + hy;endx = [a:hx:b]; y = [c:hy:d];if nargout == 4, if nargin == 8        error('Exact solution not available');    else        error = max(max(abs(u-ue)))/max(max(abs(ue)));end, endreturn