数值计算B大作业Word格式.docx
《数值计算B大作业Word格式.docx》由会员分享，可在线阅读，更多相关《数值计算B大作业Word格式.docx（57页珍藏版）》请在冰豆网上搜索。
24.4
21.70
25.2
28.5
24.8
25.4
二、数学原理
假设有n+1个不同的节点及函数在节点上的值（x
，y
），……（x
），插值多项式有如下形式：
（1）
其中系数
（i=0,1,2……n）为特定系数，可由插值样条
（i=0,1,2……n）确定。
根据均差的定义，把x看成[a,b]上的一点,可得
f（x）=f（
）+f[
]（
）
f[x,
]=f[
]+f[x,
]（
……
，…x
]=f[x,
]+f[x,
]（x-x
综合以上式子，把后一式代入前一式，可得到：
f（x）=f[
]+f[
）+f[
）（
）+
…+f[x,
）…（x-x
）+f[x,
]
=N
（x）+
其中
N
（x）=f[
）
（2）
=f（x）-N
（x）=f[x,
（3）
=（
Newton插值的系数
（i=0,1,2……n）可以用差商表示。
一般有
[
]（k=0,1,2，……，n）（4）
把（4）代入
（1）得到满足插值条件N
（i=0,1,2，……n）的n次Newton插值多项式
（x）=f（
）+……+f[
）…（
）.
其中插值余项为：
介于
之间。
三、程序设计
function[y,A,C,L]=newdscg（X,Y,x,M）
%y为对应x的值，A为差商表，C为多项式系数，L为多项式
%X为给定节点，Y为节点值，x为待求节点
n=length（X）;
m=length（x）;
%n为X的长度
fort=1:
m
z=x（t）;
A=zeros（n,n）;
A（:
1）=Y'
;
s=0.0;
p=1.0;
q1=1.0;
c1=1.0;
forj=2:
n
fori=j:
A（i,j）=（A（i,j-1）-A（i-1,j-1））/（X（i）-X（i-j+1））;
end
q1=abs（q1*（z-X（j-1）））;
c1=c1*j;
C=A（n,n）;
q1=abs（q1*（z-X（n）））;
fork=（n-1）:
-1:
C=conv（C,poly（X（k）））;
d=length（C）;
C（d）=C（d）+A（k,k）;
end
y（k）=polyval（C,z）;
%输出y值
L（k,:
）=poly2sym（C）;
%输出多项式
>
symsM,X=[1,3,5,7];
Y=[22.5,24.4,25.2,24.8];
x=10;
[y,A,C,L]=newdscg（X,Y,x,M）
y=
21.7313
A=
22.5000000
24.40000.950000
25.20000.4000-0.13750
24.8000-0.2000-0.1500-0.0021
C=
-0.0021-0.11871.452121.1688
L=
-x^3/480-（19*x^2）/160+（697*x）/480+3387/160
4、结果分析和讨论
对于不超过三次的插值多项式，x如果选取1,3,5,7这三个点能够得到较好的三次插值多项式L=-0.0021x^3-0.1187x^2+1.4521x+21.1688。
当x=10时，也即9点30分时的温度为21.7317度，结果分析知此值应是偏小的。
对于选取不同的插值节点，能够得到不同的插值多项式，误差也不尽相同。
5、完成题目的体会与收获
牛顿插值法的重要一点就是对插值节点的选取，通过本题的编程很好的加深了对其概念的理解以及提高了应用牛顿插值法的能力，学会了运用Matlab软件对牛顿插值法相关问题进行编程求解，对Matlab计算方法与程序编辑更加熟悉。
使我对这类问题的理解有了很大的提升。
题目二：
曲线拟合
在某钢铁厂冶炼过程中，根据统计数据的含碳量与时间关系，试用最小二乘法拟合含碳量与时间t的拟合曲线，并绘制曲线拟合图形。
t（分）
0510152025303540455055
01.272.162.863.443.874.154.374.514.584.024.64
2、数学原理
从整体上考虑近似函数
同所给数据点
（i=0,1,…,m）误差
（i=0,1,…,m）的大小，常用的方法有以下三种：
一是误差
（i=0,1,…,m）绝对值的最大值
，即误差向量
的∞—范数；
二是误差绝对值的和
，即误差向量r的1—范数；
三是误差平方和
的算术平方根，即误差向量r的2—范数；
前两种方法简单、自然，但不便于微分运算，后一种方法相当于考虑2—范数的平方，因此在曲线拟合中常采用误差平方和
来度量误差
（i=0，1，…，m）的整体大小。
数据拟合的具体作法是：
对给定数据
（i=0,1,…，m），在取定的函数类
中,求
使误差
（i=0,1,…,m）的平方和最小，即
=
从几何意义上讲，就是寻求与给定点
（i=0,1,…,m）的距离平方和为最小的曲线
。
函数
称为拟合函数或最小二乘解，求拟合函数
的方法称为曲线拟合的最小二乘法。
在曲线拟合中，函数类
可有不同的选取方法。
x=[0,5,10,15,20,25,30,35,40,45,50,55];
y=[0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.02,4.64]*10^（-4）;
p=polyfit（x,y,2）
p=
1.0e-04*
-0.00240.20370.2305
plot（x,y,'
r'
通过最小二乘法得到的曲线拟合多项式是p=（-0.0024x^2+0.2037x+0.2305）*10-4。
利用Matlab软件调用最小二乘法函数即可得到多项式拟合函数，由于数据较少得到的拟合曲线不够光滑，可能会存在一定的误差。
拟合曲线总体呈现增函数趋势，与数据较为吻合。
曲线拟合较常用的方法之一就包括最小二乘法，因此能够熟练使用Matlab进行数据的曲线拟合变得尤为重要。
通过完成本题的拟合过程，对于最小二乘法曲线拟合的操作更加的熟练，能够较好的完成各类数据的拟合。
题目三：
线性方程组求解
分别利用LU分解；
平方根法与改进平方根法；
追赶法求解下列几个不同类型的线性方程组，并与准确值比较结果，分析误差产生原因。
（1）设线性方程组
将A分解为一个下三角矩阵L和一个上三角矩阵U，即：
A=LU,
其中L=
，U=
解Ax=b的问题就等价于要求解两个三角形方程组：
⑴Ly=b,求y;
⑵Ux=y,求x。
设A为非奇异矩阵,且有分解式A=LU,L为单位下三角阵,U为上三角阵。
L,U的元素可以有n步直接计算定出。
用直接三角分解法解Ax=b（要求A的所有顺序主子式都不为零）的计算公式：
①
i=2,3,…，n.
计算U的第r行，L的第r列元素（i=2,3,…,n）：
②
，i=r,r+1,…,n;
③
i=r+1,…,n,且r≠n.
求解Ly=b，Ux=y的计算公式；
④
⑤
functionf1=LUresolve（a,b）;
[n,n]=size（a）;
%行数
z=size（b）%b的行数
l=[];
%l矩阵
u=[];
%u矩阵
forj=1:
1:
u（1,j）=a（1,j）;
fori=2:
l（i,1）=a（i,1）/u（1,1）;
（n-1）
sum1=0;
fork=1:
（i-1）
sum1=l（i,k）*u（k,i）+sum1;
u（i,i）=a（i,i）-sum1;
sum2=0;
sum3=0;
forj=（i+1）:
sum2=sum2+l（i,k）*u（k,j）;
sum3=sum3+l（j,k）*u（k,i）;
u（i,j）=a（i,j）-sum2;
l（j,i）=（a（j,i）-sum3）/u（i,i）;
fori=1:
l（i,i）=1;
l（i,n）=0;
l（n,n）=1;
sum4=0;
fork=1:
sum4=sum4+l（n,k）*u（k,n）;
u（n,n）=a（n,n）-sum4;
l%输出l矩阵
u%输出u矩阵
y=l\b%输出向量y
x=u\y%输出向量x
a=[42-3-1210000
86-5-3650100
42-2-132-1031
0-215-13-1194
-426-167-3323
86-8571726-35
02-13-425301
1610-11-917342-122
462-713920124
00-18-3-24-863-1];
b=[5;
12;
3;
2;
46;
13;
38;
19;
-21];
LUresolve（a,b）
z=
101
l=
1.0000000000000
2.00001.000000000000
1.000001.00000000000
0-2.00003.00001.0000000000
-1.00001.00004.0000-0.46671.000000000
2.00001.0000-5.0000-0.2667-1.64411.00000000
0-1.00003.0000-0.0667-0.08470.29691.0000000
4.0000-1.00006.0000-0.2667-3.76272.73033.18631.0000
升级会员 