用列主元消元法解线性方程组ax=b的matlab程序。
function[ra,rb,n,x]=liezhu(a,b)
b=[a b];n=length(b);ra=rank(a);
rb=rank(b);zhica=rb-ra;
if zhica>0,disp('请注意：因为ra~=rb，所以此方程组无解。')
return
endif ra==rb
if ra==n
disp('请注意：因为ra=rb=n，所以此方程组有唯一解。')
x=zeros(n,1); c=zeros(1,n+1);
for p= 1:n-1
y,j]=max(abs(b(p:n,p)))c=b(p,:)
b(p,:)b(j+p-1,:)b(j+p-1,:)c;
for k=p+1:n
m=b(k,p)/b(p,p);
b(k,p:n+1)=b(k,p:n+1)-m*b(p,p:n+1);
endend
b=b(1:n,n+1);a=b(1:n,1:n);x(n)=b(n)/a(n,n);
for q=n-1:-1:1
x(q)=(b(q)-sum(a(q,q+1:n)*x(q+1:n)))a(q,q);
endelse
disp('请注意：因为ra=rb end
end习题3.3
3.在matlab工作窗口输入程序。
> a=[1 1 1;1 3 9;1 7 49];b=[6;5;2];[ra,rb,n,x]=liezhu(a,b)
运行后输出结果。
请注意：因为ra=rb=n，所以此方程组有唯一解。
ra = 3
rb = 3
n = 3x = 6.3750
将矩阵a进行直接lu分解的matlab程序。
function hl=zhjlu(a)
n n]=size(a);ra=rank(a);
if ra~=n
disp('请注意：因为a的n阶行列式hl等于零，所以a不能进行lu分解。a的秩ra如下：')ra,hl=det(a);
return
endif ra==n
for p=1:n
h(p)=det(a(1:p,1:p));
endhl=h(1:n);
for i=1:n
if h(1,i)==0
disp('请注意：因为a的r阶主子式等于零，所以a不能进行lu分解。a的秩ra和各阶顺序主子式值hl依次如下：')hl;ra
return
endend
if h(1,i)~=0
disp('请注意：因为a的各阶主子式都不等于零，所以a能进行lu分解。a的秩ra和各阶顺序主子式值hl依次如下：')
for j=1:n
u(1,j)=a(1,j);
endfor k=2:n
for i=2:n
for j=2:n
l(1,1)=1;l(i,i)=1;
if i>j
l(1,1)=1;l(2,1)=a(2,1)/u(1,1);l(i,1)=a(i,1)/u(1,1);
l(i,k)=(a(i,k)-l(i,1:k-1)*u(1:k-1,k))/u(k,k);
elseu(k,j)=a(k,j)-l(k,1:k-1)*u(1:k-1,j);
endend
endend
hl;ra,u,l
endend
习题3.41.(1) 在matlab工作窗口输入程序。
> a=[2 4 -6;1 5 3;1 3 2];hl=zhjlu(a)
运行后输出结果。
请注意：因为a的各阶主子式都不等于零，所以a能进行lu分解。a的秩ra和各阶顺序主子式值hl依次如下：
ra = 3
u =2.0000 4.0000 -6.0000
l =1.0000 0 0
hl =2 6 18
2) 在matlab工作窗口输入程序。
> a=[1 1 6;-1 2 9;1 -2 3];hl=zhjlu(a)
运行后输出结果。
请注意：因为a的各阶主子式都不等于零，所以a能进行lu分解。a的秩ra和各阶顺序主子式值hl依次如下：
ra =3u = 1.0000 1.0000 6.0000
l = 1.0000 0 0
hl = 1 3 36
用p范数讨论ax=b解和a的性态的matlab程序。
function acp=zpjwc(a,ja,b,jb,p)
acp=cond(a,p);da=det(a);x=a\b;
dertaa=a-ja;
pnda=norm(dertaa,p);dertab=b-jb;
pndb=norm(dertab,p);
if pndb>0
jx=a\jb;pnb=norm(b,p);pnjx=norm(jx,p);dertax=x-jx;
pnjdx=norm(dertax,p);jxx=pnjdx/pnjx;pnjx=norm(jx,p);
pnx=norm(x,p);jxx=pnjdx/pnjx;xx=pnjdx/pnx;
pndb=norm(dertab,p);
xab=pndb/pnb;pnbj=norm(jb,p);xabj=pndb/pnbj;
xgxx=acp*xab;
endif pnda>0
jx=ja\b;dertax=x-jx;pnx=norm(x,p);
pnjdx=norm(dertax,p);
pnjx=norm(jx,p);jxx=pnjdx/pnjx;xx=pnjdx/pnx;
pnja=norm(ja,p);pna=norm(a,p);
pnda=norm(dertaa,p);xabj=pnda/pnja;xab=pnda/pna;
xgxx=acp*xab;
endif (acp>50)&&da<0.1)
disp('请注意：ax=b是病态的，a的p条件数acp,a的行列式值da，解x，近似解jx，解的相对误差jxx，解的相对误差估计值xgxx，b或a的相对误差xab依次如下：')
acp,da,x,jx',xx',jxx',xgxx',xab',xabj'
elsedisp('请注意： ax=b是良态的，a的p条件数acp,a的行列式值da，解x，近似解jx，解的相对误差jxx，解的相对误差估计值xgxx，b或a的相对误差xab依次如下：')
acp,da,x',jx',xx',jxx',xgxx',xab',xabj'
end习题3.6
1.在matlab工作窗口输入程序。
>a=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10];ja=a;b=[32 23 33 31]';jb=[32.1 22.9 22.
2 30.9]';acp=zpjwc(a,ja,b,jb,1)
运行后输出结果。
请注意：ax=b是良态的，a的p条件数acp,a的行列式值da，解x，近似解jx，解的相对误差jxx，解的相对误差估计值xgxx，b或a的相对误差xab依次如下：
acp =4.4880e+03da =
ans =ans =
xx =jxx =
xgxx =
xab =
xabj =
acp =4.4880e+03
用雅可比迭代解线性方程组ax=b的matlab主程序。
function x=jacdd(a,b,x0,p,wucha,max1)
n m]=size(a);
for j=1:m
a(j)=sum(abs(a(:,j)))2*(abs(a(j,j)))
endfor i=1:n
if a(i)>=0
disp('请注意：系数矩阵a不是严格对角占优的，此雅可比迭代不一定收敛')
return
endend
if a(i)<0
disp('请注意：系数矩阵a是严格对角占优的，此方程组有唯一解，且雅可比迭代收敛')
endfor k=1:max1
k;for j=1:m
x(j)=(b(j)-a(j,[1:j-1,j+1:m])*x0([1: j-1,j+1:m]))a(j,j);
endx;djwcx=norm(x'-x0,p); xdwcx=djwcx/(norm(x',p)+eps); x0=x';x1=a\b;
if (djwcx disp('请注意：雅可比迭代收敛，此方程组的精确解jx和近似解x如下：')
breakend
endif (djwcx>wucha)&&xdwcx>wucha)
数值分析第二章
第二章插值法。1 当时，求的二次插值多项式。解 则二次拉格朗日插值多项式为。2 给出的数值表。用线性插值及二次插值计算的近似值。解 由 知，若采用线性插值法计算即，则。若采用二次插值法计算时，3 给全的函数表，步长若函数表具有5位有效数字，研究用线性插值求近似值时的总误差界。解 求解近似值时，误差可...
数值分析第二章
列选主元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 triangul...
数值分析第二章作业
1.当，时，且。求的二次牛顿插值多项式 证明 2.当时，且求的牛顿插值多项式，并给出其差商型余项表达式。3.已知函数满足，求其三次牛顿插值多项式，并证明存在使得。4.已知，求的四次牛顿插值多项式，并给出其导数 假设函数的五阶导数存在 型余项表达式。5.设，为互异实数，是以其为节点的lagrange插...