数值代数上机实验报告
数值分析上机实验报告3
实验报告三题目：函数逼近——曲线拟合目的：掌握曲线拟合基本使用方法数学原理：[P,S]=polyfit(x,y,3)其中x，y为取样值，3为得出的结果的最高次数。
P为对应次数的系数，S为误差值向量，其中x,y是等长的向量，P是一个长度为m+1的向量。
结果分析和讨论：23.观察物体的运动，得出时间t与距离s的关系如表，求运动方程。
t=[0,0.9,1.9,3.0,3.9,5.0];s=[0,10,30,50,80,110];[P,S]=polyfit(t,s,5)P =-0.5432 6.4647 -26.5609 46.1436 -13.2601 -0.0000S =R: [6x6 double]df: 0normr: 1.2579e-012所以得到方程为：-13.2601x46.1436x-26.5609x6.4647x-0.5432x2345++24.在某化学反应堆里，根据实验所得分解物的质量分数y与时间t的关系，用最小拟合求y=F（t）；>> x=0:5:55;y=[0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.62,4.64];>> [P,S]=polyfit(x,y,5)P =0.0000 -0.0000 0.0002 -0.0084 0.2851 0.0082S =R: [6x6 double]df: 6normr: 0.0487所以得到方程为：0082.02851.00084.00002.023++-xxx结论：在23题中计算的结果误差为4.5769，而在24中计算的结果误差为0.0487，说明对于曲线拟合来说，总会有误差，因为取样点并不是都过拟合的曲线的。
数值分析2024上机实验报告
数值分析2024上机实验报告数值分析是计算数学的一个重要分支，它研究如何用数值方法来解决数学问题。
在数值分析的学习过程中，学生需要通过上机实验来巩固理论知识，并学会使用相应的数值方法来解决实际问题。
本篇报告将详细介绍2024年度数值分析上机实验的内容和结果。
一、实验内容2024年度数值分析上机实验分为四个部分，分别是：方程求根、插值与拟合、数值积分和常微分方程的数值解。
1.方程求根这部分实验要求使用数值方法求解给定的非线性方程的根。
常见的数值方法有二分法、牛顿法、割线法等。
在实验过程中，我们需要熟悉这些数值方法的原理和实现步骤，并对不同方法的收敛性进行分析和比较。
2.插值与拟合这部分实验要求使用插值和拟合方法对给定的一组数据进行拟合。
插值方法包括拉格朗日插值、牛顿插值等；拟合方法包括最小二乘拟合、多项式拟合等。
在实验中，我们需要熟悉插值和拟合方法的原理和实现步骤，并对不同方法的精度和稳定性进行比较。
3.数值积分这部分实验要求使用数值方法计算给定函数的积分。
常见的数值积分方法有梯形法则、辛普森法则、龙贝格积分等。
在实验过程中，我们需要熟悉这些数值积分方法的原理和实现步骤，并对不同方法的精度和效率进行比较。
4.常微分方程的数值解这部分实验要求使用数值方法求解给定的常微分方程初值问题。
常见的数值方法有欧拉法、改进的欧拉法、四阶龙格-库塔法等。
在实验中，我们需要熟悉这些数值解方法的原理和实现步骤，并对不同方法的精度和稳定性进行比较。
二、实验结果在完成2024年度数值分析上机实验后，我们得到了以下实验结果：1.方程求根我们实现了二分法、牛顿法和割线法，并对比了它们的收敛速度和稳定性。
结果表明，割线法的收敛速度最快，但在一些情况下可能会出现振荡；二分法和牛顿法的收敛速度相对较慢，但稳定性较好。
2.插值与拟合我们实现了拉格朗日插值和最小二乘拟合，并对比了它们的拟合效果和精度。
结果表明，拉格朗日插值在小区间上拟合效果较好，但在大区间上可能出现振荡；最小二乘拟合在整体上拟合效果较好，但可能出现过拟合。
数值计算方法上机实验报告
数值计算方法上机实验报告
一、实验目的
本次实验的主要目的是熟悉和掌握数值计算方法，学习梯度下降法的
原理和实际应用，熟悉Python语言的编程基础知识，掌握Python语言的
基本语法。
二、设计思路
本次实验主要使用的python语言，利用python下的numpy，matplotlib这两个工具，来实现数值计算和可视化的任务。
1. 首先了解numpy的基本使用方法，学习numpy的矩阵操作，以及numpy提供的常见算法，如矩阵分解、特征值分解等。
2. 在了解numpy的基本操作后，可以学习matplotlib库中的可视化
技术，掌握如何将生成的数据以图表的形式展示出来。
3. 接下来就是要学习梯度下降法，首先了解梯度下降法的主要原理，以及具体的实际应用，用python实现梯度下降法给出的算法框架，最终
可以达到所期望的优化结果。
三、实验步骤
1. 熟悉Python语言的基本语法。
首先是熟悉Python语言的基本语法，学习如何使用Python实现变量
定义，控制语句，函数定义，类使用，以及面向对象编程的基本概念。
2. 学习numpy库的使用方法。
其次是学习numpy库的使用方法，学习如何使用numpy库构建矩阵，学习numpy库的向量，矩阵操作，以及numpy库提供的常见算法，如矩阵分解，特征值分解等。
3. 学习matplotlib库的使用方法。
数值分析第一次上机练习实验报告
数值分析第一次上机练习实验报告一、实验目的本次实验旨在通过上机练习，加深对数值分析方法的理解，并掌握实际应用中的数值计算方法。
二、实验内容1. 数值计算的基本概念和方法在本次实验中，我们首先回顾了数值计算的基本概念和方法。
数值计算是一种通过计算机进行数值近似的方法，其包括近似解的计算、误差分析和稳定性分析等内容。
2. 方程求解的数值方法接下来，我们学习了方程求解的数值方法。
方程求解是数值分析中非常重要的一部分，其目的是找到方程的实数或复数解。
我们学习了二分法、牛顿法和割线法等常用的数值求解方法，并对它们的原理和步骤进行了理论学习。
3. 插值和拟合插值和拟合是数值分析中常用的数值逼近方法。
在本次实验中，我们学习了插值和拟合的基本原理，并介绍了常见的插值方法，如拉格朗日插值和牛顿插值。
我们还学习了最小二乘拟合方法，如线性拟合和多项式拟合方法。
4. 数值积分和数值微分数值积分和数值微分是数值分析中的两个重要内容。
在本次实验中，我们学习了数值积分和数值微分的基本原理，并介绍了常用的数值积分方法，如梯形法和辛卜生公式。
我们还学习了数值微分的数值方法，如差商法和牛顿插值法。
5. 常微分方程的数值解法常微分方程是物理和工程问题中常见的数学模型，在本次实验中，我们学习了常微分方程的数值解法，包括欧拉法和四阶龙格-库塔法。
我们学习了这些方法的步骤和原理，并通过具体的实例进行了演示。
三、实验结果及分析通过本次实验，我们深入理解了数值分析的基本原理和方法。
我们通过实际操作，掌握了方程求解、插值和拟合、数值积分和数值微分以及常微分方程的数值解法等数值计算方法。
实验结果表明，在使用数值计算方法时，我们要注意误差的控制和结果的稳定性。
根据实验结果，我们可以对计算结果进行误差分析，并选择适当的数值方法和参数来提高计算的精度和稳定性。
此外，在实际应用中，我们还需要根据具体问题的特点和条件选择合适的数值方法和算法。
四、实验总结通过本次实验，我们对数值分析的基本原理和方法有了更加深入的了解。
数值分析上机实验报告
数值分析上机实验报告摘要：本报告是对数值分析课程上机实验的总结和分析，涵盖了多种算法和数据处理方法，通过对实验结果的分析，探究了数值计算的一般过程和计算的稳定性。
1. 引言数值计算是数学的一个重要分支，广泛应用于物理、金融、工程等领域。
本次实验是对数值分析课程知识的实际应用，通过上机实现算法，探究数值计算的可靠性和误差分析。
2. 实验方法本次实验中，我们实现了多种算法，包括：（1）牛顿迭代法求方程的根；（2）高斯消元法求线性方程组的解；（3）最小二乘法拟合数据点；（4）拉格朗日插值法估计函数值；（5）梯形公式和辛普森公式求积分近似值。
对于每个算法，我们都进行了多组数值和不同参数的实验，并记录了相关数据和误差。
在实验过程中，我们着重考虑了算法的可靠性和计算的稳定性。
3. 实验结果与分析在实验中，我们得到了大量的实验数据和误差分析，通过对数据的展示和分析，我们得到了以下结论：（1）牛顿迭代法求解非线性方程的根能够对算法的初始值和迭代次数进行适当的调整，从而达到更高的稳定性和可靠性。
（2）高斯消元法求解线性方程组的解需要注意到矩阵的奇异性和精度的影响，从而保证计算的准确性。
（3）最小二乘法拟合数据点需要考虑到拟合的函数形式和数据的误差范围，采取适当的数据预处理和拟合函数的选择能够提高计算的准确性。
（4）拉格朗日插值法估计函数值需要考虑到插值点的选择和插值函数的阶数，防止出现龙格现象和插值误差过大的情况。
（5）梯形公式和辛普森公式求积分近似值需要考虑到采样密度和拟合函数的选择，从而保证计算的稳定性和收敛速度。
4. 结论通过本次实验的分析和总结，我们得到了深入的认识和理解数值计算的一般过程和算法的稳定性和可靠性，对于以后的数值计算应用也提供了一定的指导和参考。
数值代数上机实验报告
试验项目名称：平方根法与改进平方根法实验内容：先用你熟悉的计算机语言将平方根法和改进平方根法编写成通用的子程序，然后用你编写的程序求解对称正定方程组Ax=b，其中，A=[101 10 1…1 10 11 10]100*100b随机生成，比较计算结果，评论方法优劣。
实验要求：平方根法与改进的平方根的解法步骤；存储单元，变量名称说明；系数矩阵与右端项的生成；结果分析。
实验报告姓名：罗胜利班级：信息与计算科学0802 学号：u200810087实验一、平方根法与改进平方根法先用你所熟悉的计算机语言将平方根法和改进的平方根法编成通用的子程序，然后用你编写的程序求解对称正定方程组AX=b，其中系数矩阵为40阶Hilbert矩阵，即系数矩阵A的第i行第j列元素为=，向量b的第i个分量为=.平方根法函数程序如下：function [x,b]=pingfanggenfa(A,b)n=size(A);n=n(1);x=A^-1*b; %矩阵求解disp('Matlab自带解即为x');for k=1:nA(k,k)=sqrt(A(k,k));A(k+1:n,k)=A(k+1:n,k)/A(k,k);for j=k+1:n;A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);endend %Cholesky分解for j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n); %前代法A=A';for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1); %回代法disp('平方根法的解即为b');endfunction [x]=ave(A,b,n) %用改进平方根法求解Ax=b L=zeros(n,n); %L为n*n矩阵D=diag(n,0); %D为n*n的主对角矩阵S=L*D;for i=1:n %L的主对角元素均为1L(i,i)=1;for i=1:nfor j=1:n %验证A是否为对称正定矩阵if (eig(A)<=0)|(A(i,j)~=A(j,i)) %A的特征值小于0或A非对称时，输出wrong disp('wrong');break;endendendD(1,1)=A(1,1); %将A分解使得A=LDL Tfor i=2:nfor j=1:i-1S(i,j)=A(i,j)-sum(S(i,1:j-1)*L(j,1:j-1)');L(i,1:i-1)=S(i,1:i-1)/D(1:i-1,1:i-1);endD(i,i)=A(i,i)-sum(S(i,1:i-1)*L(i,1:i-1)');endy=zeros(n,1); % x，y为n*1阶矩阵x=zeros(n,1);for i=1:ny(i)=(b(i)-sum(L(i,1:i-1)*D(1:i-1,1:i-1)*y(1:i-1)))/D(i,i); %通过LDy=b解得y的值endfor i=n:-1:1x(i)=y(i)-sum(L(i+1:n,i)'*x(i+1:n)); %通过L T x=y解得x的值改进平方根法函数程序如下：function b=gaijinpinfanggenfa(A,b)n=size(A);n=n(1);v=zeros(n,1);for j=1:nfor i=1:j-1v(i)=A(j,i)*A(i,i);endA(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1);A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1))/A(j,j);end %LDL'分解B=diag(A);D=zeros(n);for i=1:nD(i,i)=B(i);A(i,i)=1;EndA=tril(A); %得到L和Dfor j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n); %前代法A=D*(A');for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1); %回代法disp('改进平方根法解得的解即为b');end调用函数解题：clear;clc;n=input('请输入矩阵维数：');b=zeros(n,1);A=zeros(n);for i=1:nfor j=1:nA(i,j)=1/(i+j-1);b(i)=b(i)+1/(i+j-1);endend %生成hilbert矩阵[x,b]=pingfanggenfa(A,b)b=gaijinpinfanggenfa(A,b)运行结果：请输入矩阵维数：40Warning: Matrix is close to singular or badly scaled.Results may be inaccurate. RCOND = 6.570692e-020.> In pingfanggenfa at 4In qiujie at 10Matlab自带解即为x平方根法的解即为bx =1.60358.96850.85621.01950.9375-50.2500-3.0000-16.000024.0000-49.5000-30.000039.000022.0000-64.0000-12.00002.000010.2500-10.5000-1.0000-10.875083.000046.0000-98.000012.0000-69.000068.000021.0000-50.7188-8.7500-8.0000 112.00006.0000 -68.750022.000044.0000 -28.00008.0000 -44.000012.0000b =1.0e+007 *0.0000-0.00000.0001-0.0004-0.00140.0424-0.29801.1419-2.73354.2539-4.30182.7733-1.19890.5406-0.36880.3285-0.44380.4621-0.25130.05650.0000-0.00510.0071-0.0027-0.0031-0.00190.00090.0002-0.0002-0.00060.00040.0001-0.00020.00010.0000-0.00000.0000-0.0000-0.0000改进平方根法解得的解即为bb =1.0e+024 *0.0000-0.00000.0001-0.00120.0139-0.09540.4208-1.21012.0624-1.0394-3.33436.2567-0.2463-7.45942.80303.69900.7277-1.7484-0.4854-3.60100.25325.18621.44100.8738-4.56541.04224.0920-2.7764-2.2148-0.89530.36654.89671.04160.1281-4.3387-1.1902-2.83348.4610-3.6008实验二、利用QR分解解线性方程组：利用QR分解解线性方程组Ax=b，其中A=[16 4 8 4;4 10 8 4;8 8 12 10;4 4 10 12];b=[32 26 38 30];求解程序如下：定义house函数：function [v,B]=house(x)n=length(x);y=norm(x,inf);x=x/y;Q=x(2:n)'*x(2:n);v(1)=1;v(2:n)=x(2:n);if n==1B=0;elsea=sqrt(x(1)^2+Q);if x(1)<=0v(1)=x(1)-a;elsev(1)=-Q/(x(1)+a);endB=2*v(1)^2/(Q+v(1)^2);v=v/v(1);endend进行QR分解：clear;clc;A=[16 4 8 4;4 10 8 4;8 8 12 10;4 4 10 12];b=[32 26 38 30];b=b';x=size(A);m=x(1);n=x(2);d=zeros(n,1);for j=1:n[v,B]=house(A(j:m,j));A(j:m,j:n)=(eye(m-j+1)-B*(v')*v)*A(j:m,j:n);d(j)=B;if j<mA(j+1:m,j)=v(2:m-j+1);endend %QR分解R=triu(A); %得到R D=A;I=eye(m,n);Q=I;for i=1:nD(i,i)=1;endH=tril(D);M=H';for i=1:nN=I-d(i)*H(1:m,i)*M(i,1:m);Q=Q*N;end %得到Qb=(Q')*b; %Q是正交阵for j=n:-1:2b(j)=b(j)/R(j,j);b(1:j-1)=b(1:j-1)-b(j)*R(1:j-1,j);endb(1)=b(1)/R(1,1); %回带法运行结果如下：R =18.7617 9.8072 15.7769 11.08640 9.9909 9.3358 7.53410 0 5.9945 9.80130 0 0 -0.5126Q =0.8528 -0.4368 -0.2297 -0.17090.2132 0.7916 -0.4594 -0.34170.4264 0.3822 0.2844 0.76890.2132 0.1911 0.8095 -0.5126b=1.000000000000001.000000000000010.9999999999999881.00000000000001实验三、Newton下山法解非线性方程组：3x-cos（yz）-=0，-81+sinz+1.06=0，exp（-xy）+20z+=0；要求满足数值解=满足或.定义所求方程组的函数：Newtonfun.mfunction F = Newtonfun(X)F(1,1)=3*X(1)-cos(X(2)*X(3))-1/2;F(2,1)=X(1)^2-81*(X(2)+0.1)^2+sin(X(3))+1.06;F(3,1)=exp(-X(1)*X(2))+20*X(3)+(10*pi-3)/3;End向量求导：Xiangliangqiudao.mfunction J=xiangliangqiudao()syms x y zX=[x,y,z];F=[3*X(1)-cos(X(2)*X(3))-1/2;X(1)^2-81*(X(2)+0.1)^2+sin(X(3))+1.06;exp(-X(1)*X(2))+20*X(3)+(10 *pi-3)/3];J=jacobian(F,[x y z]);End代值函数：Jacobi.mfunction F=Jacobi(x)F=[ 3,x(3)*sin(x(2)*x(3)), x(2)*sin(x(2)*x(3));2*x(1), -162*x(2)-81/5,cos(x(3));-x(2)/exp(x(1)*x(2)),-x(1)/exp(x(1)*x(2)),20];End方程组求解：format long; %数据表示为双精度型X1=[0,0,0]';eps=10^(-8);k=1;i=1;X2=X1-Jacobi(X1)^(-1)*Newtonfun(X1);while (norm(subs(X2-X1,pi,3.1415926),2)>=eps)&&(norm(Newtonfun(X1),2)>=eps) if norm(Newtonfun(X2),2)<norm(Newtonfun(X1),2) %判断先后两次迭代的大小X1=X2;B=inv(Jacobi(X2));C=Newtonfun(X2);X2=X2-B*C;i=i+1;elsev=1/(2^k); %引入下山因子X1=X2;B=inv(Jacobi(X2));C=Newtonfun(X2);X2=X2-v*B*C;k=k+1;endendj=i+k-1 %迭代次数X=X2 %输出结果运行结果如下：j =5X =0.500000000000000 -0.000000000000000 -0.523598775598299。
数值代数(第08周)实验报告()
A=sparse(i,j,s,n,n,3*n-2);
spy(A);
x0=A\b;
[x]=Zuiganfa(a,c,d,b);
ifnorm(x'-x0) > 1.0e-6
error('误差太大');
end
2.2.3实验（或测试）结果
得分：
??? Error using ==> Untitled2 at 14
(2)计算矩阵 的克拉特(Crout)分解
其中
【注：在实际计算中，将数组alpha和beta分别保存在数组b和c中.】
(3)解
(4)解
2.2程序清单
2.2.1写出追赶法求解三对角方程组的函数代码
得分：
function[x]=Zuiganfa(a,c,d,b);
n=length(a);
n1=length(c);
数值代数（第08周）实验报告
学号
20150
姓名
实验目的
追赶法实验
1、实验要求
给定一个 级主对角占优的三对角矩阵 ，和一个 的列向量 ，用追赶法求解线性方程组 ，其中
注：尝试用规模尽量大（即n尽量大）的矩阵来完成实验.
2、追赶法
2.1算法描述
得分：
（写出追赶法求解三对角方程组的计算过程）
(1)将矩阵 的三条对角线和常数向量 分别记为数组
%U(i,i+1)=q(i);
%end
end
2.2.2实验（或测试）代码
得分：
n=100000;
c=rand(1,n)*100;
d=rand(1,n)*100;
a = rand(1,n)*1000+c+d;
数值分析上机实验报告
数值分析上机实验报告导言：本次上机实验主要是针对数值分析课程中的一些基本算法进行实验验证。
实验内容包括迭代法、插值法、数值积分和常微分方程的数值解等。
在实验过程中，我们将会使用MATLAB进行算法的实现，并对结果进行分析。
一、迭代法迭代法是解决函数零点、方程解等问题的常用方法。
我们将选择几个常见的函数进行迭代求根的实验。
（1）二分法二分法是一种简单而有效的迭代求根法。
通过函数在区间两个端点处的函数值异号来确定函数在区间内存在零点，并通过不断缩小区间来逼近零点。
（2）牛顿法牛顿法利用函数的一阶导数和二阶导数的信息来逼近零点。
通过不断迭代更新逼近值，可以较快地求得零点。
实验结果表明，对于简单的函数，这两种迭代法都具有很好的收敛性和稳定性。
但对于一些复杂的函数，可能会出现迭代失效或者收敛速度很慢的情况。
二、插值法插值法是在给定一些离散数据点的情况下，通过构造一个插值函数来逼近未知函数的值。
本实验我们将使用拉格朗日插值和牛顿插值两种方法进行实验。
（1）拉格朗日插值拉格朗日插值通过构造一个多项式函数来逼近未知函数的值。
该多项式经过离散数据点，并且是唯一的。
该方法简单易懂，但插值点越多，多项式次数越高，插值函数的精度也就越高。
（2）牛顿插值牛顿插值利用差商的概念，通过构造一个插值多项式来逼近未知函数的值。
与拉格朗日插值相比，牛顿插值的计算过程更加高效。
但同样要求插值点的选择要合理，否则可能出现插值函数不收敛的情况。
实验结果表明，这两种插值方法都能够很好地逼近未知函数的值。
插值点的选择对插值结果有很大的影响，过多或者过少的插值点都可能导致插值结果偏离真实函数的值。
三、数值积分数值积分是一种将定积分问题转化为数值求和的方法。
本实验我们将使用复合梯形求积法和复合辛普森求积法进行实验。
（1）复合梯形求积法复合梯形求积法将定积分区间等分为若干小区间，然后使用梯形公式对每个小区间进行近似求积，最后将结果相加得到整个定积分的近似值。
数值代数(第05周)实验报告()
2.1算法描述
得分：
（写出列主元Gauss消去法的计算过程）
(1)找出矩阵A第i到m行主元绝对值最大的位置并将最大元素所在行与该列第一行互换
(2)单位矩阵做同样的行变换
(3)构成L的下三角矩阵非对角线元素做同样的行变换
2.2程序清单
2.2.1写出列主元Gauss消去法的函数代码
数值代数（第05周）实验报告
学号
201505
姓名
实验目的
学习列主元Gauss消去法
1、实验要求
1.1给定一个n x n级矩阵A，用列主元Gauss消去法将矩阵分解为PA=LU形式，其中P为一个n x n级置换矩阵，L为一个n x n级单位下三角形矩阵，U为一个n x n级上三角形矩阵.
1.2综合运用列主元Gauss消去法求解线性方程组Ax=b.
P=eye(n);
L=eye(n);
L1=zeros(n);
U=doubleபைடு நூலகம்A);
fork=1:n-1
forr=k+1:n
A(r,k)=A(r,k)/A(k,k);
A(r,k+1:n)= A(r,k+1:n)+A(r,k)*A(r,k+1:n);
end
end
U=triu(A);
L=tril(A,-1)+eye(n);
得分：
function [L,U,P] = GaussElimColPiv(A)
%列主元的Gauss消去法
%输入
% A: nxn级矩阵
%输出
% L: nxn级单位下三角矩阵形矩阵
% U: nxn级上三角形矩阵
% P:置换矩阵
数值线性代数第二版徐树方高立张平文上机习题第二章实验报告
（1）估计5到20阶Hilbert 矩阵的∞数条件数（2）设n n R A ⨯∈⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡------=111111111011001，先随机地选取n R x ∈，并计算出x A b n =；然后再用列主元Gauss 消去法求解该方程组，假定计算解为∧x 。
试对n 从5到30估计计算解∧x 的精度，并且与真实相对误差作比较。
解（1）分析：利用for 使n 从5循环到20，利用()hilb 函数得到Hilbert 矩阵A ；先将算法2.5.1编制成通用的子程序，利用算法2.5.1编成的子程序)(B opt v =，对TA B -=求解，得到∞-1A的一个估计值v v =~；再利用inf),(A norm 得到∞A ；则条件数inf),(1A norm v A A K *==∞∞-。
另，矩阵A 的∞数条件数可由inf),(A cond 直接算出，两者可进行比较。
程序为1 算法2.5.1编成的子程序)(B opt v =function v=opt(B)k=1;n=length(B);x=1./n*ones(n,1); while k==1 w=B*x;v=sign(w); z=B'*v;if norm(z,inf)<=z'*x v=norm(w,1); k=0; elsex=zeros(n,1);[s,t]=max(abs(z)); x(t)=1; k=1; end end end2 问题（1）求解 ex2_1for n=5:20A=hilb(n);B=inv(A.');v=opt(B);K1=v*norm(A,inf);K2=cond(A,inf);disp(['n=',num2str(n)])disp(['估计条件数为',num2str(K1)])disp(['实际条件数为',num2str(K2)])end计算结果为n=5估计条件数为943656实际条件数为943656n=6估计条件数为29070279.0028实际条件数为29070279.0028n=7估计条件数为985194887.5079实际条件数为985194887.5079n=8估计条件数为.7717实际条件数为.7717n=9估计条件数为86.422实际条件数为86.422n=10估计条件数为750.67实际条件数为750.67n=11估计条件数为49344实际条件数为49344Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.547634e-17.> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.547634e-17.> In cond at 47In ex2_1 at 6n=12估计条件数为3.3713e+16实际条件数为3.3713e+16Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.847381e-19.Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.847381e-19.> In cond at 47In ex2_1 at 6n=13估计条件数为1.5327e+18实际条件数为1.5327e+18Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.246123e-18.> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.246123e-18.> In cond at 47In ex2_1 at 6n=14估计条件数为4.8374e+17实际条件数为4.8374e+17Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.491876e-19.> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.491876e-19.> In cond at 47In ex2_1 at 6n=15估计条件数为4.9674e+17实际条件数为5.3619e+17Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 9.137489e-19.> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 9.137489e-19.> In cond at 47In ex2_1 at 6n=16估计条件数为8.3166e+17实际条件数为8.3167e+17Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.244518e-19.> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.244518e-19.> In cond at 47n=17估计条件数为1.093e+18 实际条件数为1.093e+18Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.693737e-19. > In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.693737e-19. > In cond at 47 In ex2_1 at 6 n=18估计条件数为2.0651e+18 实际条件数为2.7893e+18Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.264685e-19. > In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.264685e-19. > In cond at 47 In ex2_1 at 6 n=19估计条件数为2.9357e+18 实际条件数为2.9357e+18Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.351364e-19. > In ex2_1 at 3Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.351364e-19. > In cond at 47 In ex2_1 at 6 n=20估计条件数为2.674e+18 实际条件数为6.473e+18结果分析随着矩阵阶数增加，估计值误差开始出现，20,17,16,15=n 时估计条件数与实际值存在误差；且条件数很大，Hilbert 矩阵为病态的。
数值计算方法上机实验报告1
数值计算方法上机实验报告1华北电力大学上机实验报告课程名称：数值计算方法专业班级：学生姓名：学号：指导教师：张建成实验目的：复习和巩固数值计算方法的基本数学模型，全面掌握运用计算机进行数值计算的具体过程及相关问题。
利用计算机语言独立编写、调试数值计算方法程序，培养学生利用计算机和所学理论知识分析解决实际问题的能力。
上机练习任务：利用计算机基本C 语言编写并调试一系列数值方法计算通用程序，并能正确计算给定题目，掌握调试技能。
掌握文件使用编程技能，如文件的各类操作，数据格式设计、通用程序运行过程中文件输入输出运行方式设计等。
一、列主元素消去法求解线性方程组 1、算法原理为避免绝对值很小的元素作为主元，在每次消元之前增加一个选主元的过程，将绝对值大的元素交换到主对角线的位置。
列主元素消元法是当变换到第k 步时，从k 列的kk a 及以下的各元素中选取绝对值最大的元素，然后通过二交换将其交换到kk a 的位置上。
2、输入输出变量ija ：为系数矩阵的各个系数K ：表示到第k 步消元 3、具体算例输入增广矩阵为： 3 1 2 -3 8 2 1 3 22 3 2 1 28解得：1x =6，2x =4，3x =2；二、LU 分解法求解线性方程组1、算法原理应用高斯消去法解n 阶线性方程Ax b =经过1n -步消去后得出一个等价的上三角形方程组()()n n A x b =，对上三角形方程组用逐步回代就可以求出解来。
这个过程也可通过矩阵分解来实现。
将非奇异阵分解成一个下三角阵L 和上三角阵U 的乘积A LU =称为对矩阵A 的三角分解，又称LU 分解。
根据LU 分解,将Ax b =分解为Ly bUx y =??=?形式,简化了求解问题。
2、输入输出变量ij a 为系数矩阵元素i b 为常数矩阵系数,i j i jl u 分别为下、上三角矩阵元素 k 表示第k 步消元 3、具体算例输入增广矩阵 3 2 3 4 39 3 -2 2 14 4 2 3 43 解得： 6 5 3三、拉格朗日插值1、算法原理设函数()y f x =在区间[a,b]上有节点01,,,,n x x x 上的函数值，构造一个次数不超过n次的代数多项式1110()n n n n p x a x a x a x a --=++++ ,使 (),0,1,,i i P x y i n == 。
数值分析上机实验报告
数值分析上机实验报告数值分析上机实验报告一、引言数值分析是一门研究利用计算机进行数值计算的学科。
通过数值分析，我们可以使用数学方法和算法来解决实际问题，例如求解方程、插值和逼近、数值积分等。
本次上机实验旨在通过编程实现数值计算方法，并应用于实际问题中。
二、实验目的本次实验的目的是掌握数值计算方法的基本原理和实现过程，加深对数值分析理论的理解，并通过实际应用提高编程能力。
三、实验内容1. 数值求解方程首先，我们使用二分法和牛顿迭代法分别求解非线性方程的根。
通过编写程序，输入方程的初始值和精度要求，计算得到方程的根，并与理论解进行对比。
2. 数值插值和逼近接下来，我们使用拉格朗日插值和最小二乘法进行数据的插值和逼近。
通过编写程序，输入给定的数据点，计算得到插值多项式和逼近多项式，并绘制出插值曲线和逼近曲线。
3. 数值积分然后，我们使用梯形法和辛普森法进行定积分的数值计算。
通过编写程序，输入被积函数和积分区间，计算得到定积分的近似值，并与解析解进行比较。
四、实验步骤1. 数值求解方程（1）使用二分法求解非线性方程的根。
根据二分法的原理，编写程序实现二分法求解方程的根。
（2）使用牛顿迭代法求解非线性方程的根。
根据牛顿迭代法的原理，编写程序实现牛顿迭代法求解方程的根。
2. 数值插值和逼近（1）使用拉格朗日插值法进行数据的插值。
根据拉格朗日插值法的原理，编写程序实现数据的插值。
（2）使用最小二乘法进行数据的逼近。
根据最小二乘法的原理，编写程序实现数据的逼近。
3. 数值积分（1）使用梯形法进行定积分的数值计算。
根据梯形法的原理，编写程序实现定积分的数值计算。
（2）使用辛普森法进行定积分的数值计算。
根据辛普森法的原理，编写程序实现定积分的数值计算。
五、实验结果与分析1. 数值求解方程通过二分法和牛顿迭代法，我们成功求解了给定非线性方程的根，并与理论解进行了对比。
结果表明，二分法和牛顿迭代法都能够较好地求解非线性方程的根，但在不同的问题中，二者的收敛速度和精度可能会有所差异。
数值代数上机报告
数值代数上机报告一、 目的意义： 把矩阵A 分解成一个下三角阵与一个上三角阵的乘积，即 A=LR ，其中L 为下三角阵，R 为上三角阵，这样原线性方程组就可以化为⎩⎨⎧==y Rx b Ly 的求解问题，方便求解。
二、 算法：1) 输入系数矩阵 A ；2) 利用公式iji k ij ij ij r l a r ∑-=-=11和ii ki i k jk ji ji r r l a l /)(11∑-=-=交错进行，计算得出矩阵L 和R ；3) 回带到⎩⎨⎧==y Rx b Ly 中得出原线性方程组的解。
三、 源程序：#include <stdio.h>#include <string.h>#include <math.h>#include <stdlib.h>#define N 100main(){int i,j,k,s,n;printf("请输入系数矩阵A 的阶数：n= ");scanf("%d",&n);floata[N][N]={0},L[N][N]={0},R[N][N]={0},sigma1,sigma2,b[N],y[N],x[N];/*为L 主对角线元素赋1*/for(i=0;i<n;i++){L[i][i]=1;}printf("请输入系数矩阵A ：\n"); /*输入系数矩阵A*/报告1 Doolittle 分解for(i=0;i<n;i++){for(j=0;j<n;j++)scanf("%f",&a[i][j]);}for(k=0;k<n;k++){for(j=k;j<n;j++) /*计算矩阵R*/{sigma1=0;for(s=0;s<=k-1;s++)sigma1+=L[k][s]*R[s][j];R[k][j]=a[k][j]-sigma1;}for(i=k;i<n;i++) /*计算矩阵L*/{sigma2=0;for(s=0;s<=k-1;s++)sigma2+=L[i][s]*R[s][k];L[i][k]=(a[i][k]-sigma2)/R[k][k];}}printf("\n A矩阵为：\n");/*输出矩阵L、R*/ for(i=0;i<n;i++){for(j=0;j<n;j++)printf("%5.1f ",a[i][j]);printf("\n");}printf("\n L矩阵为：\n");for(i=0;i<n;i++){for(j=0;j<n;j++)printf("%5.1f ",L[i][j]);printf("\n");}printf("\n R矩阵为：\n");for(i=0;i<n;i++){for(j=0;j<n;j++)printf("%5.1f ",R[i][j]);printf("\n");}printf("请输入b矩阵\n",i+1);for(i=0;i<n;i++)scanf("%f",&b[i]);for(i=0;i<n;i++)/*回代法求解方程组Ly=b*/{sigma1=0;for(k=0;k<=i-1;k++)sigma1+=L[i][k]*y[k];y[i]=b[i]-sigma1;}for(i=n-1;i>=0;i--){sigma2=0;for(k=i+1;k<n;k++)sigma2+=R[i][k]*x[k];x[i]=(y[i]-sigma2)/R[i][i];}printf("解得x为：\n");for(i=0;i<n;i++)printf("%5.1f ",x[i]);printf("\n");}四、计算结果与分析：分析：运行结果与预想的结果相近，误差对结果的影响不是很大，比较理想五、 参考文献：[1]刑志栋. 矩阵数值分析. 陕西： 陕西科学技术出版社， 2005[2]谭浩强. C 语言程序设计. 北京：清华大学出版社，2005一、 目的意义： 对称正定矩阵是实践中经常遇到的一种特殊矩阵类型矩阵，由于矩阵本身的流量好兴致，使得cholesky 分解在存储和运算量上较一般消去法节省一半左右，且解的精度高，cholesky 分解方法是目前计算机求解该类问题最有效的方法之一。
数值代数上机实验报告
数值代数课程设计实验报告姓名： 班级： 学号： 实验日期：一、实验名称 代数的数值解法 二、实验环境 MATLAB7.0实验一、平方根法与改进平方根法一、实验要求：用熟悉的计算机语言将不选主元和列主元Gasuss 消元法编写成通用的子程序，然后用编写的程序求解下列方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--⨯14151515157681681681681681612321n n n n n x x x x x x 用所编的程序分别求解40、84、120阶方程组的解。
二、算法描述及实验步骤GAuss 程序如下：（1）求A 的三角分解：LU A =;（2）求解b y =L 得y ； （3）求解y x =U 得x ；列主元Gasuss 消元法程序如下： 1求A 的列主元分解：LU PA =;2求解b y P L =得y ； 3求解y x =U 得x ；三、调试过程及实验结果：%----------------方程系数---------------->> A1=Sanduijiaozhen(8,6,1,40); >> A2=Sanduijiaozhen(8,6,1,84); >> A3=Sanduijiaozhen(8,6,1,120); >> b1(1)=7;b2(1)=7;b3(1)=7;>> for i=2:39b1(i)=15;end>> b1(40)=14;>> for i=2:83b2(i)=15;end>> b2(40)=14;>> for i=2:119b1(i)=15;end>> b3(120)=14;%----------------方程解---------------->> x11=GAuss(A1,b1')>> x12=GAuss Zhu(A1,b1')>> x21=GAuss(A2,b2')>> x22=GAuss Zhu(A3,b3')>> x31=GAuss(A3,b3')>> x32=GAuss Zhu(A3,b3')运行结果：（n=40）GAuss消元法的解即为x11 =1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000列主元GAuss消元法的解即为x12 =1 1 1 1 1 1 1 1 1 1 111111111111111111111111111111六、源程序：function A=Sanduijiaozhen(a,b,c,n)%生成n阶以a,b,c为元素的三对角阵A=diag(b*ones(1,n),0)+diag(c*ones(1,n-1),1)+diag(a*ones(1,n-1),-1);function x=GAuss(A,b)n=length(b);x=b;%-------分解---------------for i=1:n-1for j=i+1:nmi=A(j,i)/A(i,i);b(j)=b(j)-mi*b(i);for k=i:nA(j,k)=A(j,k)-mi*A(i,k);endAB=[A,b];endend%-----------回代------------------x(n)=b(n)/A(n,n);for i=n-1:-1:1s=0;for j=i+1:ns=s+A(i,j)*x(j);endx(i)=(b(i)-s)/A(i,i);endfunction x=GAussZhu(A,b)n=length(b);x=b;%----------------------选主元---------------------for k=1:n-1a_max=0;for i=k:nif abs(A(i,k))>a_maxa_max=abs(A(i,k));r=i;endendif r>kfor j=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;end%--------------消元-----------------for i=k+1:nm=A(i,k)/A(k,k);for j=k:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);endendif abs(A(n,n))==0return;endAbZhu=[A,b];%----------------回代-----------------------x(n)=b(n)/A(n,n);for i=n-1:-1:1for j=i+1:nb(i)=b(i)-A(i,j)*x(j);endx(i)=b(i)/A(i,i);end实验二、平方根法与改进平方根法一、实验要求：用计算机语言将平方根法和改进的平方根法编成通用的子程序，然后用编写的程序求解对称正定方程组100阶方程组AX=b，二、算法描述及实验步骤：平方根法函数程序如下：1、求A 的Cholesky 分解：L L A T=;2、求解b y =L 得y ；3、求解y x =TL 得x ； 改进平方根法函数程序如下：1、求A 的Cholesky 分解：T=LDL A ; 2、求解b y =L 得y ； 3、求解y x =TDL 得x ；三、调试过程及实验结果：clear;clc;%----------------方程系数---------------->> A=Sanduijiaozhen(1,10,1,100); >> b(1)=11; >> for i=2:99 b(i)=12; end>> b(100)=11;>> x1=Cholesky(A,b') >> x2=GJCholesky(A,b')运行结果：平方根法的解即为 x1 =1.0000 1.00001.0000 1.00001.0000 1.00001.0000 1.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000改进平方根法解得的解即为x2 =1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000 1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00000.99991.00090.99081.09080.1010四、源程序：function x=Cholesky(A,b)n=size(A);n=n(1);% x=A^-1*b;% disp('Matlab自带解即为x');%-----------------Cholesky分解-------------------for k=1:nA(k,k)=sqrt(A(k,k));A(k+1:n,k)=A(k+1:n,k)/A(k,k);for j=k+1:n;A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);endend%------------------前代法求解Ly=b----------------------------for j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n);%-----------------回代法求解L'x=y-----------------------------A=A';for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1);disp('平方根法的解即为');function b=GJCholesky(A,b)n=size(A);n=n(1);v=zeros(n,1);%----------------------LDL'分解-----------------------------for j=1:nfor i=1:j-1v(i)=A(j,i)*A(i,i);endA(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1);A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1))/A(j,j);endB=diag(A);D=zeros(n);for i=1:nD(i,i)=B(i);A(i,i)=1;end%-------------------前代法---------------------------A=tril(A); %得到L和Dfor j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n);%-----------------回代法-----------------------------A=D*(A');for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1);disp('改进平方根法解得的解即为');实验三、二次多项式拟合一、实验要求：用计算机语言编制利用QR分解求解线性最小二乘问题的通用子程序，用编写的程序求解一个二次多项式使在残向量的范数最小的意义下拟合下面的数据t-1 -0.75 -0.5 0 0.25 0.5 0.75iy 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125i二、算法描述及实验步骤：QR分解求解程序如下：1、求A 的QR 分解;2、计算b c 11T =Q ;3、求解上三角方程1c x =R 得x ；三、调试过程及实验结果：>> t=[-1 -0.75 -0.5 0 0.25 0.5 0.75];>> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> plot(t,y,'r*');>> legend('实验数据(ti,yi)'); >> xlabel('t'), ylabel('y');>> title('二次多项式拟合的数据点(ti,yi)的散点图');运行后屏幕显示数据的散点图（略）.（3）编写下列MATLAB 程序计算)(x f 在),(i i y x 处的函数值，即输入程序 >> syms a b c>> t=[-1 -0.75 -0.5 0 0.25 0.5 0.75]; >> fi=a.*t.^2+ b.*t+c%运行后屏幕显示关于 ,,a b c 的线性方程组fi =[a-b+c,9/16*a-3/4*b+c,1/4*a-1/2*b+c,c,1/16*a+1/4*b+c,1/4*a+1/2*b+c,9/16*a+3/4*b +c]编写构造残向量2范数的MATLAB 程序>> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> fy=fi-y; fy2=fy.^2; J=sum(fy.^2); 运行后屏幕显示误差平方和如下 J=(a-b+c-1)^2+(9/16*a-3/4*b+c-13/16)^2+(1/4*a-1/2*b+c-3/4)^2+(c-1)^2+(1/16*a+1/4*b+c-21/16)^2+(1/4*a+1/2*b+c-7/4)^2+(9/16*a+3/4*b+c-37/16)^2为求,,a b c 使J 达到最小，只需利用极值的必要条件0J a ∂=∂，0J b ∂=∂，0J c∂=∂，得到关于,,a b c 的线性方程组，这可以由下面的MATLAB 程序完成，即输入程序 >> Ja1=diff(J,a); Ja2=diff(J,b); Ja3=diff(J,c);>> Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3) 运行后屏幕显示J 分别对,,a b c 的偏导数如下 Ja11 =451/128*a-63/32*b+43/8*c-887/128 Ja21 =-63/32*a+43/8*b-3/2*c-61/32Ja31 =43/8*a-3/2*b+14*c-143/8解线性方程组112131000Ja Ja Ja ===，，，输入下列程序 >> A=[451/128, -63/32, -3/2 ;-63/32,43/8,-3/2;43/8,-3/2,14]; >> B=[887/128,61/32,143/8]; >> C=B/A, f=poly2sym(C)运行后屏幕显示拟合函数f 及其系数C 如下 C =0.3081 0.8587 1.4018 f =924/2999*x^2+10301/11996*x+4204/2999 故所求的拟合曲线为2()0.30810.8581 1.4018f x x x =++四、源程序：>> t=[-1 -0.75 -0.5 0 0.25 0.5 0.75];>> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> plot(t,y,'r*');>> legend('实验数据(ti,yi)'); >> xlabel('t'), ylabel('y');>> title('二次多项式拟合的数据点(ti,yi)的散点图'); >> syms a b c>> t=[-1 -0.75 -0.5 0 0.25 0.5 0.75]; >> fi=a.*t.^2+ b.*t+c fi =[ a-b+c, 9/16*a-3/4*b+c, 1/4*a-1/2*b+c, c, 1/16*a+1/4*b+c, 1/4*a+1/2*b+c, 9/16*a+3/4*b+c]>> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> fy=fi-y; fy2=fy.^2; J=sum(fy.^2) J =(a-b+c-1)^2+(9/16*a-3/4*b+c-13/16)^2+(1/4*a-1/2*b+c-3/4)^2+(c-1)^2+(1/16*a+1/4*b+c-21/16)^2+(1/4*a+1/2*b+c-7/4)^2+(9/16*a+3/4*b+c-37/16)^2>> Ja1=diff(J,a); Ja2=diff(J,b); Ja3=diff(J,c);>> Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3)Ja11 =451/128*a-63/32*b+43/8*c-887/128Ja21 =-63/32*a+43/8*b-3/2*c-61/32Ja31 =43/8*a-3/2*b+14*c-143/8>> A=[451/128, -63/32, -3/2 ;-63/32,43/8,-3/2;43/8,-3/2,14]; >> B=[887/128,61/32,143/8];>> C=B/A, f=poly2sym(C)C =0.3081 0.8587 1.4018f =924/2999*x^2+10301/11996*x+4204/2999>>。
数值代数实验报告(1)绝对经典
数值代数实验报告Numerical Linear Algebra And ItsApplications学生所在学院：理学院学生所在班级：计算数学10-1学生姓名：戈东潮指导教师：于春肖教务处2012年12月实验一实验名称： Poisson 方程边值问题的五点差分格式 实验时间: 2012年12月13日 星期四 实验成绩： 一、实验目的：通过上机利用Matlab 数学软件实现Poisson 方程的边值问题的五点差分格式的线性代数方程组来认真解读Poisson 方程边值问题的具体思想与方法，使我们掌握得更加深刻，并且做到学习与使用并存，增加学习的实际动手性，不再让学习局限于书本和纸上，而是利用计算机学习来增加我们的学习兴趣。
二、实验内容：利用Poisson 方程来解下列问题：⎪⎩⎪⎨⎧∂∈=∈=∂∂+∂∂-ΩΩ),(,),(),(,sin sin )(y x y x u y x y x y u x u 0222222πππ 其中}{的边界是ΩΩ∂<<∈Ω,1,0),(y x y x 。
边值问题的解释是y x y x u ππsin sin ),(=（1）用例题中模型问题的方法取N=5,10（也可以取N=20）列出五点差分格式的线性代数方程组三、实验过程：系数矩阵A 的Matlab 实现函数如下：%文件名：poisson.m function T=poisson(N) for i=1:N a(i)=4; end b=diag(a); for j=1:N-1b(j,j+1)=-1;endfor j=2:Nb(j,j-1)=-1;endfor i=1:NI(i)=-1;endI=diag(I);for i=1:N:N*Nfor j=1:NT(i+j-1,i:i+N-1)=b(j,:);endendfor i=1:N:N*N-Nfor j=1:NT(i+j-1,i+N:i+N-1+N)=I(j,:);endendfor i=1+N:N:N*Nfor j=1:NT(i+j-1,i-N:i-1)=I(j,:);endend求解常数项b的Matlab实现函数如下：%函数名：constant.mfunction b=constant(N)n=N+1;h=1/n;i=1;for y=1/n:1/n:N/nfor x=1/n:1/n:N/nf(i)=2*pi*pi*sin(pi*x)*sin(pi*y);i=i+1;endendb=h*h*f;b=b';四、实验结果(总结/方案)在Matlab运行窗口输入>>A=poisson(5) Enter输出结果A =Columns 1 through 114 -1 0 0 0 -1 0 0 0 0 0-1 4 -1 0 0 0 -1 0 0 0 00 -1 4 -1 0 0 0 -1 0 0 00 0 -1 4 -1 0 0 0 -1 0 00 0 0 -1 4 0 0 0 0 -1 0-1 0 0 0 0 4 -1 0 0 0 -10 -1 0 0 0 -1 4 -1 0 0 00 0 -1 0 0 0 -1 4 -1 0 00 0 0 -1 0 0 0 -1 4 -1 00 0 0 0 -1 0 0 0 -1 4 00 0 0 0 0 -1 0 0 0 0 40 0 0 0 0 0 -1 0 0 0 -10 0 0 0 0 0 0 -1 0 0 00 0 0 0 0 0 0 0 -1 0 00 0 0 0 0 0 0 0 0 -1 00 0 0 0 0 0 0 0 0 0 -10 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 Columns 12 through 220 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 00 -1 0 0 0 0 0 0 0 0 00 0 -1 0 0 0 0 0 0 0 00 0 0 -1 0 0 0 0 0 0 0 -1 0 0 0 -1 0 0 0 0 0 04 -1 0 0 0 -1 0 0 0 0 0 -1 4 -1 0 0 0 -1 0 0 0 00 -1 4 -1 0 0 0 -1 0 0 00 0 -1 4 0 0 0 0 -1 0 00 0 0 0 4 -1 0 0 0 -1 0 -1 0 0 0 -1 4 -1 0 0 0 -10 -1 0 0 0 -1 4 -1 0 0 00 0 -1 0 0 0 -1 4 -1 0 00 0 0 -1 0 0 0 -1 4 0 00 0 0 0 -1 0 0 0 0 4 -10 0 0 0 0 -1 0 0 0 -1 40 0 0 0 0 0 -1 0 0 0 -10 0 0 0 0 0 0 -1 0 0 00 0 0 0 0 0 0 0 -1 0 0 Columns 23 through 250 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 -1 0 0 0 -1 0 0 0 -1 0 0 4 -1 0 -1 4 -10 -1 4 在运行窗口输入 >>b=constant(5) Enter 常数项的输出结果b =[ 0.1371 0.2374 0.2742 0.2374 0.1371 0.2374 0.4112 0.4749 0.4112 0.2374 0.2742 0.4749 0.5483 0.4749 0.2742 0.2374 0.4112 0.4749 0.4112 0.2374 0.1371 0.2374 0.2742 0.2374 0.1371]T所以A*u=b 的五点差分格式为j i j i j i j i j i j i f h u u u u u ,,,,,,211114=-----+-+ （i,j=1,2……）实验二实验名称： 用Jacobi 迭代法和SOR 迭代法求解方程组 实验时间: 2012年12月13日 星期四 实验成绩： 一、实验目的：通过上机利用Matlab 数学软件采用Jacobi 迭代法和SOR 迭代法来求解方程组，求解过程中了解两种迭代方法的基本思想与迭代过程，并分析两种迭代方法的收敛性和实用用以及他们的异同点。
数值代数实验报告1
暨南大学本科实验报告专用纸课程名称数值代数成绩评定实验项目名称(填写实验所属章节名称) 指导教师刘娟实验项目编号0701******* 实验项目类型验证实验地点机房学生姓名王伟文学号2010051727学院信息科学技术学院数学系信息与计算科学专业2010级实验时间2011年9月 1 日～12月30日温度24℃第一章1.实验选题：（写出你在该项目所选的实验题目）第一章上机练习12.谈谈你对该算法的理解：（简单谈一下你是如何理解该算法的？）对算法的理解：先将84阶的矩阵A分解为一个下三角矩阵L和上三角矩阵U,先考虑下三角形方程组Ly=b,利用前代法解出y,再考虑上三角形方程组Ux=y;利用回代法解出x。
列主元gauss消去法Ax b PA LU Ly Pb Ux y=⇔===,,;3.实验内容（将实验程序及其实验结果粘贴，最好对程序各部分注释清楚，比如设置了哪些函数，这些函数的输入输出是什么，具有什么功能？）实验程序对矩阵A进行LU分解的程序function [ L,U ] = LUfac( A )for k=1:n-1A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);endL=tril(A,0);for i=1:nL(i,i)=1;endU=triu(A,0);end利用前代法解出y值的程序function [ b ] = TSL( L,b )n=size(L,1);for j=1:n-1b(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);endb(n)=b(n)/L(n,n);end利用回代法解出x值的程序function [ b ] = TSU( U,b )n=size(U,1);for j=n:-1:2b(j)=b(j)/U(j,j);b(1:j-1)=b(1:j-1)-b(j)*U(1:j-1,j);endb(1)=b(1)/U(1,1);end主函数程序（‘生成84阶的矩阵A’）A=eye(84);A=6*A;for i=2:84A(i,i-1)=8;A(i-1,i)=1;end（‘生成84乘1的矩阵b’）b=ones(84,1);b=b*15;b(1)=7;b(84)=14;[L,U]=LUfac(A);（‘调用函数LUfac对矩阵A进行分解’）y=TSL(L,b);（‘调用函数TSL求解Ly=b方程’）x=TSU(U,y);（‘调用函数TSU求解UX=b方程’）用MATLAB求解得到的结果x’ans =1.0e+008 *Columns 1 through 70.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 8 through 140.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 15 through 210.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 22 through 280.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 29 through 350.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 36 through 420.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 43 through 490.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 50 through 560.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 Columns 57 through 63-0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 Columns 64 through 700.0000 -0.0000 0.0000 -0.0001 0.0002 -0.0003 0.0007 Columns 71 through 77-0.0013 0.0026 -0.0052 0.0105 -0.0209 0.0419 -0.0836 Columns 78 through 840.1665 -0.3303 0.6501 -1.2582 2.3487 -4.0263 5.3684列主元gauss消去法function [L,U,P]=Lufac(A)n=size(A,1);P=eye(n,n);for i=1:n-1[r,m]=max(abs(A(i:n,i)));m=m+i-1;A([i,m],:)=A([m,i],:);P([i,m])=P([m,i]);if A(i,i)~=0A(i+1:n,i)=A(i+1:n,i)/ A(i,i);A(i+1:n,i+1:n)= A(i+1:n,i+1:n)-A(i+1:n,i)*A(i,i+1:n); endendU=triu(A);L=tril(A,-1)+eye(n,n);endA=eye(84);A=6*A;for i=2:84A(i,i-1)=8; A(i-1,i)=1; endb=ones(84,1);b=b*15;b(1)=7;b(84)=14;[L,U,P]=Lufac(A); b=P*b;y=TSL(L,b);x=TSU(U,y);求解结果为x =1.0e+025 *0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00010.0002-0.00040.0007-0.00150.0030-0.00600.0119-0.02380.0475-0.09460.1878-0.36960.7154-1.33552.2894-3.05251.教师评语、评分：（请认真对待每次数值代数实验，期末交实验报告，将计算实验成绩）2.实验选题：第三章上机习题1（3）。
数值代数实验报告
数值代数实验报告数值代数实验报告引言：数值代数是一门研究数值计算方法和算法的学科，它在科学计算和工程应用中起着重要的作用。
本实验报告旨在通过实际的数值计算问题，探讨数值代数的应用和效果。
实验一：线性方程组求解线性方程组求解是数值代数中的一个重要问题。
在实验中，我们使用了高斯消元法和LU分解法两种求解线性方程组的方法，并对比了它们的效果。
首先，我们考虑一个3×3的线性方程组：2x + 3y - z = 54x - 2y + 2z = 1x + y + z = 3通过高斯消元法，我们将该方程组转化为上三角形式，并得到解x=1, y=2, z=0。
而通过LU分解法，我们将该方程组分解为LU两个矩阵的乘积，并得到相同的解。
接下来，我们考虑一个更大的线性方程组，例如10×10的方程组。
通过比较高斯消元法和LU分解法的运行时间，我们可以发现LU分解法在处理大规模方程组时更加高效。
实验二：特征值与特征向量计算特征值与特征向量计算是数值代数中的另一个重要问题。
在实验中，我们使用了幂法和QR方法两种求解特征值与特征向量的方法，并对比了它们的效果。
首先，我们考虑一个3×3的矩阵：1 2 34 5 67 8 9通过幂法，我们可以得到该矩阵的最大特征值为15.372，对应的特征向量为[0.384, 0.707, 0.577]。
而通过QR方法，我们也可以得到相同的结果。
接下来，我们考虑一个更大的矩阵，例如10×10的矩阵。
通过比较幂法和QR 方法的运行时间，我们可以发现QR方法在处理大规模矩阵时更加高效。
实验三：奇异值分解奇异值分解是数值代数中的一种重要技术，它可以将一个矩阵分解为三个矩阵的乘积，从而实现数据降维和信息提取的目的。
在实验中，我们使用了奇异值分解方法，并通过实际的数据集进行了验证。
我们选取了一个包含1000个样本和20个特征的数据集，通过奇异值分解，我们将该数据集分解为三个矩阵U、S和V的乘积。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值代数课程设计实验报告姓名： 班级： 学号： 实验日期：一、实验名称 代数的数值解法 二、实验环境 MATLAB7.0实验一、平方根法与改进平方根法一、实验要求：用熟悉的计算机语言将不选主元和列主元Gasuss 消元法编写成通用的子程序，然后用编写的程序求解下列方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--⨯14151515157681681681681681612321n n n n n x x x x x x 用所编的程序分别求解40、84、120阶方程组的解。
二、算法描述及实验步骤 GAuss 程序如下：（1）求A 的三角分解：LU A =;（2）求解b y =L 得y ； （3）求解y x =U 得x ；列主元Gasuss 消元法程序如下： 1求A 的列主元分解：LU PA =;2求解b y P L =得y ； 3求解y x =U 得x ；三、调试过程及实验结果：%----------------方程系数---------------->> A1=Sanduijiaozhen(8,6,1,40); >> A2=Sanduijiaozhen(8,6,1,84); >> A3=Sanduijiaozhen(8,6,1,120); >> b1(1)=7;b2(1)=7;b3(1)=7;>> for i=2:39b1(i)=15;end>> b1(40)=14;>> for i=2:83b2(i)=15;end>> b2(40)=14;>> for i=2:119b1(i)=15;end>> b3(120)=14;%----------------方程解---------------->> x11=GAuss(A1,b1')>> x12=GAuss Zhu(A1,b1')>> x21=GAuss(A2,b2')>> x22=GAuss Zhu(A3,b3')>> x31=GAuss(A3,b3')>> x32=GAuss Zhu(A3,b3')运行结果：（n=40）GAuss消元法的解即为x11 =1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000列主元GAuss消元法的解即为x12 =1 1 1 1 1 1 1 1 1 1 111111111111111111111111111111六、源程序：function A=Sanduijiaozhen(a,b,c,n)%生成n阶以a,b,c为元素的三对角阵A=diag(b*ones(1,n),0)+diag(c*ones(1,n-1),1)+diag(a*ones(1,n-1),-1);function x=GAuss(A,b)n=length(b);x=b;%-------分解---------------for i=1:n-1for j=i+1:nmi=A(j,i)/A(i,i);b(j)=b(j)-mi*b(i);for k=i:nA(j,k)=A(j,k)-mi*A(i,k);endAB=[A,b];endend%-----------回代------------------x(n)=b(n)/A(n,n);for i=n-1:-1:1s=0;for j=i+1:ns=s+A(i,j)*x(j);endx(i)=(b(i)-s)/A(i,i);endfunction x=GAussZhu(A,b)n=length(b);x=b;%----------------------选主元---------------------for k=1:n-1a_max=0;for i=k:nif abs(A(i,k))>a_maxa_max=abs(A(i,k));r=i;endendif r>kfor j=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;end%--------------消元-----------------for i=k+1:nm=A(i,k)/A(k,k);for j=k:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);endendif abs(A(n,n))==0return;endAbZhu=[A,b];%----------------回代-----------------------x(n)=b(n)/A(n,n);for i=n-1:-1:1for j=i+1:nb(i)=b(i)-A(i,j)*x(j);endx(i)=b(i)/A(i,i);end实验二、平方根法与改进平方根法一、实验要求：用计算机语言将平方根法和改进的平方根法编成通用的子程序，然后用编写的程序求解对称正定方程组100阶方程组AX=b，二、算法描述及实验步骤：平方根法函数程序如下：1、求A 的Cholesky 分解：L L A T=; 2、求解b y =L 得y ；3、求解y x =TL 得x ；改进平方根法函数程序如下：1、求A 的Cholesky 分解：T=LDL A ; 2、求解b y =L 得y ； 3、求解y x =TDL 得x ；三、调试过程及实验结果：clear;clc;%----------------方程系数---------------->> A=Sanduijiaozhen(1,10,1,100); >> b(1)=11; >> for i=2:99 b(i)=12; end>> b(100)=11; >> x1=Cholesky(A,b') >> x2=GJCholesky(A,b')运行结果：平方根法的解即为 x1 =1.0000 1.00001.0000 1.00001.0000 1.00001.0000 1.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000改进平方根法解得的解即为x2 =1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000 1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00000.99991.00090.99081.09080.1010四、源程序：function x=Cholesky(A,b)n=size(A);n=n(1);% x=A^-1*b;% disp('Matlab自带解即为x');%-----------------Cholesky分解-------------------for k=1:nA(k,k)=sqrt(A(k,k));A(k+1:n,k)=A(k+1:n,k)/A(k,k);for j=k+1:n;A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);endend%------------------前代法求解Ly=b----------------------------for j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n);%-----------------回代法求解L'x=y-----------------------------A=A';for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1);disp('平方根法的解即为');function b=GJCholesky(A,b)n=size(A);n=n(1);v=zeros(n,1);%----------------------LDL'分解-----------------------------for j=1:nfor i=1:j-1v(i)=A(j,i)*A(i,i);endA(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1);A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1))/A(j,j);endB=diag(A);D=zeros(n);for i=1:nD(i,i)=B(i);A(i,i)=1;end%-------------------前代法---------------------------A=tril(A); %得到L和Dfor j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n);%-----------------回代法-----------------------------A=D*(A');for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1);disp('改进平方根法解得的解即为');实验三、二次多项式拟合一、实验要求：用计算机语言编制利用QR分解求解线性最小二乘问题的通用子程序，用编写的程序求解一个二次多项式使在残向量的范数最小的意义下拟合下面的数据t-1 -0.75 -0.5 0 0.25 0.5 0.75iy 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125i二、算法描述及实验步骤：QR分解求解程序如下：1、求A 的QR 分解;2、计算b c 11T=Q ;3、求解上三角方程1c x =R 得x ；三、调试过程及实验结果：>> t=[-1 -0.75 -0.5 0 0.25 0.5 0.75];>> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> plot(t,y,'r*');>> legend('实验数据(ti,yi)'); >> xlabel('t'), ylabel('y');>> title('二次多项式拟合的数据点(ti,yi)的散点图');运行后屏幕显示数据的散点图（略）.（3）编写下列MATLAB 程序计算)(x f 在),(i i y x 处的函数值，即输入程序 >> syms a b c>> t=[-1 -0.75 -0.5 0 0.25 0.5 0.75]; >> fi=a.*t.^2+ b.*t+c%运行后屏幕显示关于 ,,a b c 的线性方程组fi =[a-b+c,9/16*a-3/4*b+c,1/4*a-1/2*b+c,c,1/16*a+1/4*b+c,1/4*a+1/2*b+c,9/16*a+3/4*b +c]编写构造残向量2范数的MATLAB 程序>> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> fy=fi-y; fy2=fy.^2; J=sum(fy.^2); 运行后屏幕显示误差平方和如下 J=(a-b+c-1)^2+(9/16*a-3/4*b+c-13/16)^2+(1/4*a-1/2*b+c-3/4)^2+(c-1)^2+(1/16*a+1/4*b+c-21/16)^2+(1/4*a+1/2*b+c-7/4)^2+(9/16*a+3/4*b+c-37/16)^2为求,,a b c 使J 达到最小，只需利用极值的必要条件0J a ∂=∂，0J b ∂=∂，0Jc∂=∂，得到关于,,a b c 的线性方程组，这可以由下面的MATLAB 程序完成，即输入程序 >> Ja1=diff(J,a); Ja2=diff(J,b); Ja3=diff(J,c);>> Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3) 运行后屏幕显示J 分别对,,a b c 的偏导数如下 Ja11 =451/128*a-63/32*b+43/8*c-887/128 Ja21 =-63/32*a+43/8*b-3/2*c-61/32Ja31 =43/8*a-3/2*b+14*c-143/8解线性方程组112131000Ja Ja Ja ===，，，输入下列程序 >> A=[451/128, -63/32, -3/2 ;-63/32,43/8,-3/2;43/8,-3/2,14]; >> B=[887/128,61/32,143/8]; >> C=B/A, f=poly2sym(C)运行后屏幕显示拟合函数f 及其系数C 如下 C =0.3081 0.8587 1.4018 f =924/2999*x^2+10301/11996*x+4204/2999 故所求的拟合曲线为2()0.30810.8581 1.4018f x x x =++四、源程序：>> t=[-1 -0.75 -0.5 0 0.25 0.5 0.75];>> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> plot(t,y,'r*');>> legend('实验数据(ti,yi)'); >> xlabel('t'), ylabel('y');>> title('二次多项式拟合的数据点(ti,yi)的散点图'); >> syms a b c>> t=[-1 -0.75 -0.5 0 0.25 0.5 0.75]; >> fi=a.*t.^2+ b.*t+c fi =[ a-b+c, 9/16*a-3/4*b+c, 1/4*a-1/2*b+c, c, 1/16*a+1/4*b+c, 1/4*a+1/2*b+c, 9/16*a+3/4*b+c]>> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> fy=fi-y; fy2=fy.^2; J=sum(fy.^2) J =(a-b+c-1)^2+(9/16*a-3/4*b+c-13/16)^2+(1/4*a-1/2*b+c-3/4)^2+(c-1)^2+(1/16*a+1/4*b+c-21/16)^2+(1/4*a+1/2*b+c-7/4)^2+(9/16*a+3/4*b+c-37/16)^2>> Ja1=diff(J,a); Ja2=diff(J,b); Ja3=diff(J,c);>> Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3)Ja11 =451/128*a-63/32*b+43/8*c-887/128Ja21 =-63/32*a+43/8*b-3/2*c-61/32Ja31 =43/8*a-3/2*b+14*c-143/8>> A=[451/128, -63/32, -3/2 ;-63/32,43/8,-3/2;43/8,-3/2,14]; >> B=[887/128,61/32,143/8];>> C=B/A, f=poly2sym(C)C =0.3081 0.8587 1.4018f =924/2999*x^2+10301/11996*x+4204/2999>>word文档可自由复制编辑。