浙大计算方法上机报告.docx
《浙大计算方法上机报告.docx》由会员分享，可在线阅读，更多相关《浙大计算方法上机报告.docx（17页珍藏版）》请在冰点文库上搜索。
浙大计算方法上机报告
学号：
**********姓名：
专业：
作业1：
用列主元高斯消去法和列主元三角分解法解P227页第3题
1.列主元高斯消去法
目的：
用高斯消去法解Ax=b时，其中设A为非奇异矩阵，可能出现a
=0情况，这时必须进行带行交换的高斯消去法。
但在实际计算中即使
但其绝对值很小时，用
作除数，会导致中间结果矩阵A（k）元素数量级严重增长和舍入误差的扩散，使得最后的计算结果不可靠。
列主元高斯消去法可以难过一般高斯法的这些缺点。
一、列主元高斯消去法解方程的Matlab程序如下：
functiona=columneli（a）%对矩阵a进行列主元消去
[mn]=size（a）;%求取a的行数m和列数n
fori=1:
m-1
[maxEle,pos]=max（abs（a（i:
end,i）））;
maxRow=pos+i-1;%在每次变换前寻找绝对值最大主所在列maxRow
ifa（maxRow,i）==0
disp（'矩阵为奇异矩阵'）
return%对于非奇异矩阵，在程序中给予提示,结束程序
end
if（maxRow~=i）
temp=a（maxRow,:
）;a（maxRow,:
）=a（i,:
）;a（i,:
）=temp;
end%与列主元绝对值最大的行进行行交换
forj=i+1:
m
a（j,i）=a（j,i）/a（i,i）;%求取第j列主元
fork=i+1:
n
a（j,k）=a（j,k）-a（j,i）*a（i,k）;%对第j列主元进行行变换
end
end
end
functionx=elisolve（a）%利用列主元消去的结果求方程的解,a为方程组的增广矩阵
a=columneli（a）;
[m,n]=size（a）;
x=zeros（m,1）;
fori=m:
-1:
1
x（i）=（a（i,n）-a（i,i:
m）*x（i:
m））/a（i,i）;
end
二、列主元高斯消去法解P227页第3题：
答案：
x=[7/6-1/31/2]T
三、程序流程图如下：
条件图框里面没有条件式，因为i是从m到1，所以i=1后运行下面的命令，直接输出结果。
这里也是一样的，直接运行下面的命令。
四、Matlab中的命令截图如下：
运行结果：
x=[7/6-1/31/2]T
上机结果和正确答案完全相同。
2.列主元三角分解法
一，列主元三角分解法解方程的Matlab程序如下：
function[l,u,p]=tridecom（a）%对矩阵a进行三角分解
[mn]=size（a）;%求取a的行数m和列数n
if（m~=n）
disp（'矩阵的行数与列数不相等，无法进行三角分解'）
return
end
l=zeros（n）;u=zeros（n）;p=eye（n）;%p*a=l*u
[maxEle,goodRow]=max（abs（a（:
1）））;
temp=a（1,:
）;a（1,:
）=a（goodRow,:
）;a（goodRow,:
）=temp;
temp=p（1,:
）;p（1,:
）=p（goodRow,:
）;p（goodRow,:
）=temp;%行交换
u（1,:
）=a（1,:
）;l（1:
end,1）=a（1:
end,1）/u（1,1）;
forr=2:
n
tempUrr=zeros（n,1）;
fork=r:
n
tempUrr（k）=a（k,r）-l（k,1:
r-1）*u（1:
r-1,r）;
end
[maxEle,goodRow]=max（abs（tempUrr））;
temp=a（r,:
）;a（r,:
）=a（goodRow,:
）;a（goodRow,:
）=temp;
temp=l（r,:
）;l（r,:
）=l（goodRow,:
）;l（goodRow,:
）=temp;
temp=p（r,:
）;p（r,:
）=p（goodRow,:
）;p（goodRow,:
）=temp;
fori=r:
n
u（r,i）=a（r,i）-l（r,1:
r-1）*u（1:
r-1,i）;
end
fori=r:
n
l（i,r）=（a（i,r）-l（i,1:
r-1）*u（1:
r-1,r））/u（r,r）;
end
end
functionx=lusolve（a）%用三角分解法解线性方程组Ax=b，a为方程组的增广矩阵[A|b]
[m,n]=size（a）;
[l,u,p]=tridecom（a（:
1:
m））;%对A进行三角分解A=p*l*u
x=p*a（:
n）;%对a的最后一列进行与方程组系数三角分解时相同的行变换
fori=2:
m
x（i）=x（i）-l（i,1:
i-1）*x（1:
i-1）;%求出（l^-1）*p*b
end
x（m）=x（m）/u（m,m）;
fori=m-1:
-1:
1
x（i）=（x（i）-u（i,i+1:
m）*x（i+1:
m））/u（i,i）;%求出方程组的解x=（U^-1）*（L^-1）*p*b
end
二、程序流程图如下：
二、解P227页第3题在Matlab中的截图如下：
可见，p*a=l*u
作业2；编写Matlab程序，用最小二乘法拟合多项式曲线
（解P89页第1题）
目的和意义；根据观测或实验得到一系列的数据，确定了与自变量的某些点相应的函数值。
当函数比较复杂或根本无法写出解析式时，往往根据观测数据构造一个适当的简单的函数近似地代替要寻求的函数。
利用最小二乘法拟合多项式来描述其函数。
一、用最小二乘法拟合多项式曲线的Matlab程序如下：
functiona=polyercheng（x,y,n）
x=x';y=y';
m=size（x,1）;
c=ones（m,1+n）;
fori=1:
n
c（:
i+1）=x.^i;
end
A=c'*c;b=c'*y;
a=（A\b）';
a=fliplr（a）;
c=polyval（a,x）-y;
d=max（abs（c））;
fprintf（'最大偏差：
%f\n',d）;
d=sqrt（c'*c）;
fprintf（'均方误差：
%f',d）;
b=linspace（min（x）,max（x）,500）;
c=polyval（a,b）;
plot（x,y,'g*',b,c,'r'）;
end
二、程序的流程图如下：
三、最小二乘法拟合多项式曲线解P89页第1题:
已知实验数据如下：
x
2
4
6
8
Y
2
11
28
40
用最小二乘法求一次和二次拟合多项式，分别算出均方误差和最大偏差。
结果：
一次拟合多项式
最大偏差：
2.700000
均方误差：
3.271085
拟合多项式曲线：
y=6.5500x-12.5000
二次拟合多项式
最大偏差：
1.950000
均方误差：
2.906888
拟合多项式曲线：
y=0.1875x2+4.6750x-8.7500
四、解P89页第1题在Matlab中的截图如下：
y=6.5500x-12.5000
y=0.1875x2+4.6750x-8.7500
其中a为拟合多项式的系数阵列。
通过最小二乘法拟合多项式得到的多项式接近与实验值（代入后）。
作业3：
用龙贝格算法计算,
使截断误差不超过0.5e-6。
目的和意义；在一元函数的积分学中,我们已经熟知,若函数f（x）在区间［a,b］上连续且其原函数为F（x）,则可用牛顿―莱布尼兹公
来求定积分。
牛顿―莱布尼兹公式虽然在理论上或在解决实际问题中都起了很大的作用,但它并不能完全解决定积分的计算问题。
当被积函数f（x）在区间［a,b］上连续时,要使得复合梯形公式比较精确地代替定积分
可将分点（即基点）加密,也就是将区间［a,b］细分,然后利用复合梯形公式求积
将收敛缓慢的梯形值序列加工成收敛迅速的
，这种加速方法称为龙贝格算法。
它可以把复杂的积分比较精确的求出来。
一、龙贝格算法Matlab程序如下：
functionz=romberg（ff,a,b,epsi,k）%用龙贝格算法计算积分
ifnargin==4
k=16;
end%如果只有4个传入参数，则将二分次数k设置为默认值16
t=zeros（1,k+1）;
formatlong;
functiony=f（x）%嵌套函数，用来判断积分区域内的奇点
y=ff（x）;%若ff参数为@（x）sin（x）/x,则有y=sin（x）/x
ifisnan（y）==true||isinf（y）==true
disp（'请确保所给函数在积分区域内无奇点'）;
end
end
t
（1）=（b-a）/2*（f（a）+f（b））;
forl=1:
k
sum=0;
fori=1:
（2^（l-1））
sum=sum+f（a+（（2*i-1）*（b-a）/（2^l）））;
end
t（l+1）=t（l）/2+（b-a）/（2^l）*sum;
end%计算梯形序列
form=1:
k
temp=t;
fori=m+1:
k+1%将新序列放入t（m+1）--t（k+1）中
t（i）=（（4^m）*temp（i）-temp（i-1））/（4^m-1）;
end
ifabs（t（m+1）-t（m））break;
end
end
z=t（m+1）;
fprintf（'截断误差为%d\n',t（m+1）-t（m））;
end
二、用龙贝格算法计算P120页第5题：
计算
使截断误差不超过0.5e-6。
结果：
截断误差为6.632355e-008
=0.946083070387222
三、程序流程图如下
三、Matlab中的截图如下：
二分次数缺省时：
指定二分次数为10次时：
积分区域内有奇点时：
使用Matlab自带函数计算积分时：
升级会员 