function chap2_fdm_elliptic_1D
% 一維橢圓方程求解(常微分方程邊值問題)
% -u'' + q(x)u = f(x), 0<x<1, 取q(x)= x, f(x)=(x-1)exp(x)
% u(0)= 1, u(1)= e; 邊界條件
% 真解為 u = exp(x)N = 20;
h = 1/N;
x_all =(0:h:1)';
x = x_all(2:end-1);% 方程組右端的處理
u0 = 1;
uN = exp(1);
b = f(x);
b(1) = b(1)+u0/h^2;
b(N-1) = b(N-1)+uN/h^2;% 方程組左端矩陣的建立
B = zeros(N,N);
for i = 1:N-1B(i,i) = -2;B(i,i+1) = 1;B(i+1,i) = 1;
end
B = -1/h^2*B(1:N-1,1:N-1);D =zeros(N-1,N-1);
q = q(x);
for i = 1:N-1D(i,i) = q(i);
end
A = B+D;%求解線性方程組(Gauss消去法求解)
u = gauss(A,b);
u_e = u_exact(x_all);figure(1)
plot(x_all,[u0;u;uN],'r*',x_all,u_e,'b');
endfunction y = u_exact(x)
y = exp(x);
endfunction y = q(x)
y = x;
endfunction y = f(x)
y =(x-1).*exp(x);
endfunction x = gauss(A,B)
%消元過程
n = size(A,2);
m = size(B,2);for i = 1:mb = B(:,i);for k=1:n-1A(k+1:n,k)=A(k+1:n,k)/A(k,k); %算子A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);b(k+1:n)=b(k+1:n)-b(k)*A(k+1:n,k); end %回代過程U=A;b=b;for j =n:-1:2b(j)=b(j)/U(j,j);b(1:j-1)=b(1:j-1)-b(j)*U(1:j-1,j);endb(1)=b(1)/U(1,1);B(:,i)=b; %這個就是線性方程組的解了
end
x = B;
end
二、二維橢圓型微分方程
matlab代碼：
function chap2_fdm_possion_solution
% 二維的possion方程求解
% -Delta u = -2*exp(x+y) 邊界(a,b)*(c,d)
% u(x,0)= exp(x)
% u(x,1)= exp(x+1)
% u(0,y)= exp(y)
% u(1,y)= exp(y+1)
global mLap M_V I hx hya = 0;
b = 1;
c = 0;
d = 1;
M = 40;
N = 30;
hx =(b-a)/M;
hy =(d-c)/N;
x =(a:hx:b)';
y = (c:hy:d)';[X Y]= meshgrid(x,y);
x_in =(a+hx:hx:b-hx);
y_in =(c+hy:hy:d-hy);[X_in Y_in]= meshgrid(x_in,y_in);%生成Dxx
d_0 = -2*ones(M-1,1);
d_1 = ones(M-1,1);
d_m1 = d_1;
d_xx = spdiags([d_m1 d_0 d_1],[-1,0,1],M-1,M-1);
d_xx = d_xx/hx^2;
I_N = speye(N-1,N-1);
D_xx = kron(I_N,d_xx);%D_yy
d_0 = -2*ones(N-1,1);
d_1 = ones(N-1,1);
d_m1 = d_1;
d_yy = spdiags([d_m1 d_0 d_1],[-1,0,1],N-1,N-1);
d_yy = d_yy/hy^2;
I_M = speye(M-1,M-1);
D_yy = kron(d_yy,I_M);%左端矩陣
mLap = -(D_xx+D_yy);%右端矩陣
F = f(X_in,Y_in);
F = F';
g_d = g(x_in,0);
g_u = g(x_in,1);
g_l = g(0,y_in);
g_r = g(1,y_in);
F(:,1) = F(:,1)+g_d'/hy^2;
F(:,N-1)= F(:,N-1)+g_u'/hy^2;
F(1,:) = F(1,:)+g_l/hx^2;
F(M-1,:) = F(M-1,:)+g_r/hx^2;%求解u
b = F(:);
u = mLap\b;
U = reshape(u,M-1,N-1);
U = U';%繪制結果
u_e = u_exact(X,Y);
u_num = u_e;
u_num(2:N,2:M)= U;
figure(1)
mesh(X,Y,u_e);
figure(2)
mesh(X,Y,u_num);end%函數
function y = f(x,y)
y = -2*exp(x+y);
endfunction u_e = u_exact(x,y)
u_e = exp(x+y);
endfunction g = g(x,y)
g = exp(x+y);
end
