列选主元lu分解法:
function [l,u,pv]=luex(a)
luex : lu factorization with partial pivoting
synopsis: [l,u,pv]=luex(a)
input: a = coefficient matrix
output: l :lower triangular matrix
u ;upper triangular matrix
pv:b→pb b=b(pv)
format short
m,n]=size(a);
if m~=n, error(' a matrix needs to be square');end
pv = 1:n)';
for i=1:n-1
% partial pivoting for column i
[pivot,p] =max(abs(a(i:n,i)))
kk = p + i -1 ;
if kk ~=i
a([i kk],:a([kk i],:
pv([i kk]) pv([kk i]);
end% end of partial pivoting for column i.
% lu factorization
pivot = a(i,i);
for k = i+1:n
a(k,i)=a(k,i)/pivot;
a(k,i+1:n)=a(k,i+1:n)-a(k,i)*a(i,i+1:n);
endend
l = eye(size(a))+tril(a,-1); extract l and u
u = triu(a);
测试luex 函数:
test luex
a = 2 4 -2 -2; 1 2 4 -3; -3 -3 8 -2; -1 1 6 -3];
b = 4; 5; 7; 7];
l,u,pv]=luex(a);
y = l\b(pv);
x = u\y;
exact solution
x1 = a\b;
err=x-x1
结果如下err = 0
从结果可以看出,列主元lu分解法效果较好。
1.解:1)a =
x,flag]=gauss(a)
结果如下:x = 9
flag = 12)a =
x,flag]=gauss(a)
结果如下:x = 0.2124
flag =1
3. 解:lu分解的matlab程序:
function [l,u]=lup(a)
lup : lu factorization
synopsis: [l,u]=lup(a)
input: a = coefficient matrix
output: l :lower triangular matrix
u ;upper triangular matrix
format short
m,n]=size(a);
if m~=n, error(' a matrix needs to be square');end
pv = 1:n)';
lu factorization
for i=1:n-1
pivot = a(i,i);
for k = i+1:n
a(k,i)=a(k,i)/pivot;
a(k,i+1:n)=a(k,i+1:n)-a(k,i)*a(i,i+1:n);
endend
l = eye(size(a))+tril(a,-1); extract l and u
u = triu(a);
本题:a=[1 0 2 0; 0 1 1 1; 2 0 -1 1; 0 0 1 1];
[l,u]=lup(a)
结果如下:l =
u =4.解。a=[5 7 9 10; 6 8 10 9; 7 10 8 7; 5 7 6 5];
b=[1 1 1 1]'
l,u]=lup(a);
y = l\b;
x = u\y
结果如下:x =
5.解:a=[3 4 5; -1 3 4; -2 3 -5;];
b=[2,-2 6]';
l,u,pv]=luex(a);
y = l\b(pv);
x = u\y
结果如下:x = 1
6.解:cholesky分解法的matlab实现:
function l = cholesky(a)
cholesky : cholesky factorization of a symmetric ,positive matrix
synopsis : l = cholesky(a)
input : a = symmetric positive matrix
output : l = upper triangular matrix,and a = l'*l
m,n] =size(a);
if m~=n, error('a must be square');end
l = zeros(n);
cholesky factorization process
for i = 1:n
for j = i:n
if j==1
s = a(i,j);
else s = a(i,j) -l(1:i-1,i)'*l(1:i-1,j);
endif j > i
l(i,j) =s/l(i,i);
elseif s<=0 ,error('l is not positive definite to working precision');end
l(i,i) =sqrt(s);
endend
end对于本题:
a=[3 2 3; 2 2 0; 3 0 12;];
b = 5 3 7]’;
l=cholesky(a)
y=l'\b;
x=l\y可以得到如下结果:l =
x =对比用a\b算出的相对精确解:x = 1.0000
可以看出,上述cholesky分解法求解方程组的**结果较为准确。
9.解:x=[2 -4 3 ]'
norm(x,1)ans =
norm(x,2)ans =
norm(x,inf)ans =
11. 解:
a=[2 1 -3 -1; 3 1 0 7; -1 2 4 -2; 1 0 -1 5;];
norm(a,1)ans =
norm(a,2)ans =
norm(a,inf)
ans =11
14.解:a=[100 99;99 98];
cond(a,inf)
ans =3.9601e+04
cond(a,2)
ans =3.9206e+04
19.解:1):因为该方程组的系数矩阵a:[ 5 2 1
严格对角占有,所以对jacobi迭代法和gauss-seidel迭代法该方程组都收敛。
jacobi迭代法计算:
a=[5 2 1;-1 4 2; 2 -3 10];
b=[-12 20 3]';
x0 = zeros ( size(b) )
delta = 1e-4;
[x,iternum,flag]=jacobi(a,b,x0,delta)
结果如下:the jacobi method converges . x =
iternum =15
从上面的计算可以看出,jacobi迭代法在迭代15次后收敛,得到满足精度要求的近似解。
gauss-seidel迭代法计算:
a=[5 2 1;-1 4 2; 2 -3 10];
数值分析第二章
第二章插值法。1 当时,求的二次插值多项式。解 则二次拉格朗日插值多项式为。2 给出的数值表。用线性插值及二次插值计算的近似值。解 由 知,若采用线性插值法计算即,则。若采用二次插值法计算时,3 给全的函数表,步长若函数表具有5位有效数字,研究用线性插值求近似值时的总误差界。解 求解近似值时,误差可...
数值分析第二章作业
1.当,时,且。求的二次牛顿插值多项式 证明 2.当时,且求的牛顿插值多项式,并给出其差商型余项表达式。3.已知函数满足,求其三次牛顿插值多项式,并证明存在使得。4.已知,求的四次牛顿插值多项式,并给出其导数 假设函数的五阶导数存在 型余项表达式。5.设,为互异实数,是以其为节点的lagrange插...
数值分析作业 第二章
习题。1.当x 1,1,2时,f x 0,3,4,求的二次插值多项式。1 用单项式基底 2 用拉格朗日插值基底 3 用牛顿基底。证明三种方法得到的多项式是相同的。解 1 假设f x 的二次插值多项式为 由于 x 1,1,2时,f x 0,3,4则有 求得 则有 2 用拉格朗日插值基底 由于 则有 拉...