西安交通大学计算方法实验报告——基于matlab

西安交通大学计算方法实验报告——基于matlab内容摘要：

86018427387904(5761533065688605*(x 900)*(x 1910)*(x1920)*(x1930))/73786976294838206464+(5016235413793213*(x1900)*(x1910)*(x1920)*(x1930)*(x1940))/1180591620717411303424(5436283353602131*(x1900)*(x1910)*(x1920)*(x1930)*(x1940)*(x1950))/37778931862957161709568 + (3912362197104435*(x1900)*(x1910)*(x1920)*(x1930)*(x 1940)*(x 1950)*(x 1960))/1208925819614629174706176 (3789694605011005*(x 1900)*(x 1910)*(x 1920)*(x 1930)*(x 1940)*(x 1950)*(x 1960)*(x 1970))/77371252455336267181195264 + (8429618423594317*(x 1900)*(x 1910)*(x 1920)*(x 1930)*(x 1940)*(x 1950)*(x 1960)*(x 1970)*(x 1980))/19807040628566084398385987584 + (7234635122020997*(x 1900)*(x 1910)*(x 1920)*(x 1930)*(x 1940)*(x 1950)*(x 1960)*(x 1970)*(x 1980)*(x 1990))/5070602400912917605986812821504 (6099905420751575*(x 1900)*(x 1910)*(x 1920)*(x 1930)*(x 1940)*(x 1950)*(x 1960)*(x 1970)*(x 1980)*(x 1990)*(x 2020))/40564819207303340847894502572032 591927/200 N0 = Newton 插值 2020 年的产量 N01 = 计算方法实验报告 20 b) 三次样条插值实验结果 三次样条插值 2020 年的产量及插值函数 f = (67507170809569*t)/21990232555520+((6226409313348779*t)/2161727821137838080417169423994368193/72057594037927936)*(t2020)^212853306597599315/2199023255552 f0 = 三次样条插值 2020 年的产量 f1 = 结果分析 由实验结果可以看出， 在对 2020 年产量计算的时候， Newton 插值出现了很明显的偏差，而三次样条插值表现的很好。 与 Newton 多项式插值相比较 三次样条 插值显然得到的结果更理想，这是由于多项式插值，当 n 较大时，插值结果可能出现不一致收敛的 现象，并且高次插值计算复杂，此外三次样条插值还具有光滑性的优点，因此插值计算时应尽量应该三次样条的方法。 计算方法实验报告 21 问题 5 五、 假定某天的气温变化记录如下表所示，试用最小二乘法找出这一天的气温变化的规律 ,试计算这一天的平均气温。 时刻 0 1 2 3 4 5 6 7 8 9 10 11 12 平均气温 15 14 14 14 14 15 16 18 20 20 23 25 28 时刻 13 14 15 16 17 18 19 20 21 22 23 24 平均气温 31 32 31 29 27 25 24 22 20 18 17 16 考虑使用二次函数、三次函数、四次函数以及指数型的函数 2()b t cC ae ,计算相应的系数，估算误差，并作图比较各种函数之间的区别。 问题分析 从整体上考虑近似函数 同所给数据点 (i=0,1,„ ,m) 误差 (i=0,1,„ ,m) (i=0,1, „ ,m) 绝对值的最大值 ，即误差 向量的∞ —范数；二是误差绝对值的和 ，即误差向量 r 的 1—范数；三是误差平方和 的算术平方根，即误差向量 r 的 2—范数；前两种方法简单、自然，但不便于微分运算 ，后一种方法相当于考虑 2—范数的平方，因此在曲线拟合中常采用误差平方和 来 度量误差 (i=0， 1，„， m)的整计算方法实验报告 22 数据拟合的具体作法是：对给定数据 (i=0,1,„， m)，在取定的函数类 中 ,求 ,使误差 (i=0,1,„ ,m)的平方和最小，即 = 从几何意义上讲，就是寻求与给定点 (i=0,1,„ ,m)的距离平方和为最。 函数 称为拟合 函数或最小二乘解，求拟合函数的方法称为曲线拟合的最小二乘法。 在曲线拟合中，函数类 可有不同的选取方法 ，如多项式拟合算法。 假设给定数据点 (i=0,1,„ ,m)， 为所有次数不超过 的多项式构成的函数类，现求一 ,使得 (1) 当拟合函数为多项式时，称为多项式拟合，满足式（ 1）的 称为最小二乘拟合多项式。 特别地，当 n=1 为 的多元函数，因此上述问题即为求 的极值 问题。 由多元函数求极值的必要条件，得 (2) 即 (3) （ 3）是关于 的线性方程组，用矩阵表示为 计算方法实验报告 23 (4) 式（ 3）或式（ 4 可以证明，方程组（ 4）的系数矩阵是一个对称正定矩阵，故存在唯一解。 从式（ 4）中解出 (k=0,1,„， n) (5) 可以证明，式（ 5）中的 满足式（ 1），即 为所求的拟合多项式。 我们把 称为最小二乘拟合多项式 由式 (2) (6) 求解 的方法主要有下面两种 ， 本算法 并分别 采用 两种方法求解。 i. 求解法方程法 如果 直接利用（ 4）进行求解解方程，解法 比较简单， 但是一方面运算量较大，并且 由于条件数的影响 可能会有较大误差。 本设计采用了高斯消去法，减少了误差。 ii. 正交化法 对生成的矩阵进行 QR 分解， 上式（ 4）可写成矩阵形式， GTGa = GTy，利用矩阵的 QR 分解，可以 G = QT (RO)，（ QR 分解算法 具体实现过程见附录中的Matlab 程序 ）， Qy = (h1h2)，因此 Ra=h1，利用解方程的回代公式可以求得向量 a。 计算方法实验报告 24 具体实现程序见附录 5. 用 2()b t cC ae 拟合时，可以考虑两边取对数 log(C) = log(a) −bt2 +2bct−c2，然后对 (c,log(C))按二次多项式拟合。 Matlab 求解程序见附录五。 实验结果 a) 法方程解法实验结果 拟合二次函数如下及相应的系数 p = Px = (1238*x^2)/13455 + (2882931353438435*x)/1125899906842624 + 986/117 估算误差 e = 拟合三次函数如下及相应的系数 p = Px=(4727727282891159*x^3)/576460752303423488+(6016*x^2)/29601 (487295332134639*x)/2251799813685248 + 13072/975 估算误差 e = 拟合四次函数如下及相应的系数 p = Px=(8277009104528411*x^4)/9223372036854775808 (7389688649119097*x^3)/144115188075855872+ (3869637877579055*x^2)/4503599627370496 (8004500866520783*x)/2251799813685248 + 94306/5655 估算误差 e = 计算方法实验报告 25 拟合指数函数如下的系数 p = b = c = a = 估算误差 e = b) 正交分解法 拟合二次函数如下及相应的系数 p1 = Px1 = (1238*x^2)/13455 + (5765862706876871*x)/2251799813685248 + 986/117 估算误差 e1 = 拟合三次函数如下及相应的系数 p2 = Px2=(2363863641445577*x^3)/288230376151711744+(6016*x^2)/29601 计算方法实验报告 26 (1949181328538565*x)/9007199254740992 + 13072/975 估算误差 e2 = 拟合四次函数如下及相应的系数 p3 = Px3 = (8277009104528155*x^4)/9223372036854775808 (3694844324559451*x^3)/72057594037927936+ (7739275755157927*x^2)/9007199254740992 (8004500866520577*x)/2251799813685248 + 94306/5655 估算误差 e3 = 拟合指数函数的系数如下 p4 = b = c = a = 估算误差 e4 = 计算方法实验报告 27 结果分析 由结果可以看出，随着拟合次数的增多误差有明显减少的趋势，指数函数拟合跟二次函数拟合曲线很接近，但相对二次函数拟合误差减少了，四次函数拟合误差最小。 计算方法实验报告 28 附录 附录一 问题 1Matlab 求解 %计算 S=1/(1^2)+1/(2^2)+1/(3^2)+...+1/(100000^2),误差小于 1*10e6 clc。 clear。 sum1=0。 a=1。 %估算截取的项数，使误差小于 1*10e6 while ((1/a1/100000)(1/(a+1)1/100001))/21e6 a=a+1。 end %限定误差 disp(39。 最大误差 :39。 ) e1=((1/(a1)1/100000)(1/(a)1/100001))/2 %截取的项数 disp(39。 截取的项数 :39。 ) b=a1 %按从小到大的顺序相加 for k=b:1:1 sum1=sum1+1/(k.^2)。 end sum1=sum1+((1/(b)1/100000)+(1/(b+1)1/100001))/2。 disp(39。 估算值 :39。 ) fprintf(39。 %\n\n39。 ,sum1) %准确值，为了证明方法的可靠性 sum=0。 for k=100000:1:1 sum=sum+1/(k.^2)。 end disp(39。 准确值 :39。 ) fprintf(39。 %\n39。 ,sum) %实际误差 disp(39。 实际误差 :39。 ) ess=abs((sumsum1))/sum 计算方法实验报告 29 附录二 问题 2Matlab 求解 %对非奇异矩阵 a 作 lu 分解 clc clear disp(39。 请输入矩阵 A39。 ) A=input(39。 请输入系数矩阵 A=:39。 ) %判断被分解矩阵的维数 [n,n]=size(A)。 %给出 U 的第一行 L 的第一列 for j=1:n U(1,j)=A(1,j)。 L(j,1)=A(j,1)/U(1,1)。 end i=1。 for k=2:n for j=k:n for m=1:(k1) %计算公式中对应的 l*u 的和 uu(m)=L(k,m)*U(m,j)。 end %u 的分解结果 U。