% 北太天元 求解線性方程組?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); % 函數(shù)句柄

u_exact = @(x) sin(pi*x); % 精確解的函數(shù)句柄

x_left = 0; x_right = 1;

n = 40;

h =?(x_right -?x_left) /n ;

x = x_left:h:x_right; % 剖分網(wǎng)格

x = x ' ; % 轉(zhuǎn)成列向量

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('數(shù)值解', '精確解')

hold off