Matlab数值算法实现.ppt

Ch10. 数值算法实现
§1. 线性方程组解法
1、三角形线性方程组解法
以上三角形线性方程组为例。
回代:
%文件
function u = uptri ( a, b )
n= size(a,1) ;
x(n)= b(n) / a(n, n) ;
for i = n-1:-1:1
s=0 ;
for j = i+1: n
s=s+a(i, j) * x( j ) ;
end
x(i) = ( b(i) –s) / a(i, i) ;
end
u = x ;
求解
2、顺序Gauss消去法
(1)消去过程:
第 k 步,计算
(2)回代过程:
%文件
function u = gauss (a, b)
n = length (b) ;
for k=1: n –1
for i = k+1 : n
mult = a ( i, k) / a (k, k) ;
for j =k +1: n
% if abs ( a( k, k) ) > 1e–6
a (i, j) = a (i, j) – mult * a(k, j) ;
% else
% disp (‘顺序Gauss消去法失败’);
% pause ;
% exit ;
% end
end
可去掉%,
判断主元是否为0
b (i) = b (i) – mult * b (k) ;
end
end
x(n)= b(n) / a(n, n) ;
for i = n-1:-1:1
s=0 ;
for j = i+1: n
s=s+a(i, j) * x( j ) ;
end
x(i) = ( b(i) –s) / a(i, i) ;
end
u = x ;
回代
例:
%
a=[6, -2, 2, 4; 12, -8, 6, 10;
3,-13, 9, 3; -6, 4, 1, -18 ];
b=[16, 26, -19, -34];
x= gauss (a, b);
disp ( ‘方程组解为:’);
x
则有:
>> main
方程组解为:
x=
3 1 -2 1
3、Jacobi迭代法
%
function y = jacobi ( a, b, x0)
D = diag ( diag (a) ) ;
U = - triu (a, 1) ;
L = - tril (a, -1) ;
B = D \ ( L+U) ;
f = D \ b ;
y = B* x0 + f ; n=1;
while norm (y –x0 ) >= –6
x0 = y ;
y = B * x0 + f ; n = n +1;
end
y
n
例: P249
>> a =[10, -1, 0; -1, 10, -2; 0, -2, 10];
>> b =[9; 7; 6];
>> jacobi ( a, b, [0; 0; 0] )
y =



n =
11
解向量
迭代次数
§2. 方程求根
1、二分法
%
function y = erfen (fun,a, b, esp )
if nargin < 4 esp =1e –4 ; end
if feval (fun, a) * feval (fun, b) < 0
n = 1 ; c = (a+ b) / 2 ;
while (b-a) > esp
if feval ( fun, a) * feval (fun, c) <0
b = c ; c = ( a+b) / 2 ;
elseif feval ( fun, c) * feval (fun, b) <0
a = c ; c = ( a+b) / 2 ;
else y = c ; esp = 10000 ;
end
n= n+1 ;
end