北太天元用中心差分求解柏松方程的例子-展示稀疏矩阵S的构造和线性方程组Su=b的

% 北太天元 求解线性方程组 S u = b , 其中S 是稀疏矩阵.

%例如求解一维Poisson方程的Dirichlet边值问题

% Poisson 方程 -u'' = pi*pi* sin(pi * x) , x \in (0,1)

% 边值 u(0) = 0; u(1) = 0;

% 精确解是 u(x) = sin(pi*x)

% 

f = @(x) pi*pi* sin(pi*x); % 函数句柄

u_exact = @(x) sin(pi*x); % 精确解的函数句柄

x_left = 0; x_right = 1;

n = 40;

h = (x_right - x_left) /n ;

x = x_left:h:x_right; % 剖分网格

x = x ' ; % 转成列向量

b = h^2* f(x) ;

%边值

u0 = 0; u1= 0;

b(1) = u0;

b(end) = u1;

/*

  1  0  0  0  0

  -1 2  -1 0  0

   0  -1  2 -1 0

   0  0  -1 2  -1

   0  0  0  0  1

*/

if(n == 4)

   % 仅仅对 n = 4 的情形

   ii = [ 1 2 2 2  3 3 3 4 4 4 5 ]; 

   jj = [ 1 1 2 3  2 3 4 3 4 5 5 ];

   vv = [ 1 -1 2 -1 -1 2 -1 -1 2 -1 1 ];

   A = sparse(ii,jj,vv);

   u = A \ b

end

%对于n 取其他值的情形

iii = [ (2:n)' (2:n)' (2:n)' ];

iii = iii' ;

ii = [ 1 ; iii(:) ; n+1]

jjj = [ (1:n-1)' , (2:n)', (3:n+1)' ]

jjj = jjj';

jj = [ 1 ; jjj(:); n+1]

vvv = [-1 2 -1];

vvv = repmat(vvv, (n-1), 1);

vvv = vvv';

vv = [ 1 ; vvv(:); 1 ] ;

A = sparse(ii,jj,vv)

u = A\b

plot(x, u, 'b-o')

hold on

if(n < 100)

   xx = x_left:(x_right-x_left)/100: x_right;

else

   xx = x;

end

plot(xx, u_exact(xx), 'r-*')

legend('数值解', '精确解')

hold off