数值分析实习题作业.docx
《数值分析实习题作业.docx》由会员分享，可在线阅读，更多相关《数值分析实习题作业.docx（12页珍藏版）》请在冰豆网上搜索。
数值分析实习题作业
MATLAB数值分析程序
第一题：
编写列主元消去法程序，计算如下题目，并比较分析两者的区别。
clc;
clear;
disp（'请输入A矩阵'）;
A=[3.016.031.99;1.274.16-1.23;0.987-4.819.34];
disp（'请输入b'）;
b=[111];
A%显示系数矩阵A
b%显示b
n=length（b）;%求n
p=1:
n;
lu=A;
y=[];%记录消去法中常数项b的中间变量
fork=1:
n
[c,i]=max（abs（lu（k:
n,k）））;%找出最大值放入C，以及最大值的位置放入i
ik=i+k-1;%求最大值对应元素下标
ifik~=k
m=p（k）;p（k）=p（ik）;p（ik）=m;%记录第几行作为主元
ck=lu（k,:
）;lu（k,:
）=lu（ik,:
）;lu（ik,:
）=ck;%换主元
end
ifk==n%循环终止条件
break;
end
lu（k+1:
n,k）=lu（k+1:
n,k）/lu（k,k）;%求K+1行对应的系数
lu（k+1:
n,k+1:
n）=lu（k+1:
n,k+1:
n）-lu（k+1:
n,k）*lu（k,k+1:
n）;%高斯消去
end
l=diag（ones（n,1））+tril（lu,-1）;%生成下三角L
u=triu（lu）;%生成上三角阵U
y
（1）=b（p
（1））;
fori=2:
n
y（i）=b（p（i））-l（i,1:
i-1）*y（1:
i-1）';%计算列主元消去步骤后的常数项
end
x（n）=y（n）/u（n,n）;
fori=n-1:
-1:
1
x（i）=（y（i）-u（i,i+1:
n）*x（i+1:
n）'）/u（i,i）;%列向量存放x
end
x=x';%求出x
x
detA=det（A）;%求行列式A的值
detA
condA=cond（A）;%求条件数
condA
第
（2）题：
clc;
clear;
disp（'请输入A矩阵'）;
A=[3.006.031.99;1.274.16-1.23;0.990-4.819.34];
disp（'请输入b'）;
b=[111];
A%显示系数矩阵A
b%显示b
n=length（b）;%求n
p=1:
n;
lu=A;
y=[];%记录消去法中常数项b的中间变量
fork=1:
n
[c,i]=max（abs（lu（k:
n,k）））;%找出最大值放入C，以及最大值的位置放入i
ik=i+k-1;%求最大值对应元素下标
ifik~=k
m=p（k）;p（k）=p（ik）;p（ik）=m;%记录第几行作为主元
ck=lu（k,:
）;lu（k,:
）=lu（ik,:
）;lu（ik,:
）=ck;%换主元
end
ifk==n%循环终止条件
break;
end
lu（k+1:
n,k）=lu（k+1:
n,k）/lu（k,k）;%求K+1行对应的系数
lu（k+1:
n,k+1:
n）=lu（k+1:
n,k+1:
n）-lu（k+1:
n,k）*lu（k,k+1:
n）;%高斯消去
end
l=diag（ones（n,1））+tril（lu,-1）;%生成下三角L
u=triu（lu）;%生成上三角阵U
y
（1）=b（p
（1））;
fori=2:
n
y（i）=b（p（i））-l（i,1:
i-1）*y（1:
i-1）';%计算列主元消去步骤后的常数项
end
x（n）=y（n）/u（n,n）;
fori=n-1:
-1:
1
x（i）=（y（i）-u（i,i+1:
n）*x（i+1:
n）'）/u（i,i）;%列向量存放x
end
x=x';%求出x
x
detA=det（A）;%求行列式A的值
detA
condA=cond（A）;%求条件数
condA
结果比较：
有些矩阵，系数的微小扰动会使解向量产生巨大差别
第二题：
对龙格函数做拉格朗日多项式插值，并画出龙格现象的函数
第一步，定义原函数
functionf=f（x）
f=1./（1+25*x.^2）;
end
第二步：
计算拉格朗日差值函数
functionlangrange=langrange（x,n）
langrange=0;
xx=linspace（-1,1,n+1）;
fori=1:
n+1
lix=1;
forj=1:
n+1
ifj~=i
lix=lix.*（（x-xx（j））./（xx（i）-xx（j）））;
end
end
langrange=f（xx（i））.*lix+langrange;
end
end
第三步：
定义画龙格函数图像
functionrunge1（n）
x=linspace（-1,1,100）;
plot（x,f（x）,x,langrange（x,n））;
end
第四步：
调用上面三个函数
subplot（3,1,1）,runge1（5）,title（'5个节点'）;
subplot（3,1,2）,runge1（10）,title（'10个节点'）;
subplot（3,1,3）,runge1（20）,title（'20个节点'）;
第三题：
计算龙贝格函数积分
function[R,k,T]=romberg（f,a,b,tol）
formatlonge;
k=0;
n=1;
h=b-a;
T=h/2*（f（a）+f（b））;
err=1;
whileerr>=tol
k=k+1;
h=h/2;
tmp=0;
fori=1:
n
tmp=tmp+f（a+（2*i-1）*h）;
end
T（k+1,1）=T（k）/2+h*tmp;
forj=1:
k
T（k+1,j+1）=T（k+1,j）+（T（k+1,j）-T（k,j））/（4^j-1）;
end
n=n*2;
err=abs（T（k+1,k+1）-T（k,k））;
end
R=T（k+1,8）;
积分子函数：
functionf=f（x）
f=@（x）1./（（sin（x））^2+1/4*（cos（x））^2）;
end
主函数：
[R,k,T]=romberg（f,0,pi/2,1e-7）
合并三个函数：
functionromberg（）
functionf=f（x）
f=@（x）1./（（sin（x））^2+1/4*（cos（x））^2）;
end
function[R,k,T]=romberg（f,a,b,tol）
formatlonge;
k=0;%µü´ú´ÎÊý
n=1;%Çø¼ä¸öÊý
h=b-a;
T=h/2*（f（a）+f（b））;
err=1;
whileerr>=tol
k=k+1;
h=h/2;
tmp=0;
fori=1:
n
tmp=tmp+f（a+（2*i-1）*h）;
end
T（k+1,1）=T（k）/2+h*tmp;
forj=1:
k
T（k+1,j+1）=T（k+1,j）+（T（k+1,j）-T（k,j））/（4^j-1）;
end
n=n*2;
err=abs（T（k+1,k+1）-T（k,k））;
end
R=T（k+1,8）;
end
[R,k,T]=romberg（f,0,pi/2,1e-7）
end
升级会员 