1、数值分析报告资料数值分析课程名称 数值分析实验报告 专业班级 应用数学111班 姓名 陈伟 学号 8 教学教师 岳雪芝 实验二 函数逼近与曲线拟合 一、问题提出 从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。 在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量与时间t的拟合曲线。 t(分)0 5 10 15 20 25 30 35 40 45 50 55 0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64 二、要求 1、用最小二乘法进行曲线拟合； 2
2、、近似解析表达式为；3、打印出拟合函数，并打印出与的误差，； 4、另外选取一个近似表达式，尝试拟合效果的比较； 5、* 绘制出曲线拟合图。 三、目的和意义 1、掌握曲线拟合的最小二乘法； 2、最小二乘法亦可用于解超定线代数方程组； 3、探索拟合函数的选择与拟合精度间的关系。四、实验步骤:输入代码：t=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; for i=1:12 A(i,1)=t(i); A(i,2)=t(i).2; A(i,3)=t(i).3;
3、end p_star=Ay; y_star=p_star(1)*t+p_star(2)*t.2+p_star(3)*t.3; s_star=0; for i=1:12 s_star=s_star+(y_star(i)-y(i)2; end plot(t,y,*,t,y_star,g);比较拟和结果，分别取二次，三次和四次多项式，再做最小二乘法，代码如下： t=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; %3次拟合 p3=polyfit(t,y,3)
4、y_new=polyval(p3,t); s_3=0; for i=1:12s_3=s_3+(y_new(i)-y(i)2; end %2次拟合 p2=polyfit(t,y,2); y_new2=polyval(p2,t); s_2=0; for i=1:12 s_2=s_2+(y_new2(i)-y(i)2; end %4次拟合 p4=polyfit(t,y,4); y_new4=polyval(p4,t); s_4=0; for i=1:12 s_4=s_4+(y_new4(i)-y(i)2; end比较4种拟合函数，结论：常数项为0的三次拟合函数在保证拟合精度的同时，保证0点的函数值仍
5、为0，是四种拟合方法中最好的一种。原图与曲线拟合图如下所示：分析: 从上面的拟合中也可以知道多项式拟合误差平方和随着拟合多项式次数的增加而逐渐减少，拟合曲线更靠近实际数据，拟合更准确。实验三 数值积分与数值微分一、基本题 选用复合梯形公式，复合Simpson公式，Romberg算法，计算 （1）二、要求 1、 编制数值积分算法的程序； 2、 分别用两种算法计算同一个积分，并比较其结果； 3、 分别取不同步长，试比较计算结果（如n = 10, 20等）； 4、 给定精度要求，试用变步长算法，确定最佳步长。 三、目的和意义 1、 深刻认识数值积分法的意义； 2、 明确数值积分精度与步长的关系； 3
6、、 根据定积分的计算方法，结合专业考虑给出一个二重积分的计算问题。 四、实验步骤：代码1，复合梯形公式function result=trapezoid_integration(f,count_node,a,b) %composite trapezoid integration,f,函数表达式,count_node插入节点数量,a,开始点,b结束点 h=(b-a)/count_node;%h步长 np1=count_node+1;%包括端点,总共的点的数量 x=a:h:b;%x自变量 c=ones(np1,1); c(1,1)=0.5;c(np1,1)=0.5;%两个端点都取一半 syms s
7、ymbol_x; fx=subs(f,symbol_x,x); result=vpa(h*fx*c,9);代码2，复合simpson公式function result=simpson_integration(f,count_node,a,b) %composite simpson integration,f,函数表达式,count_node插入节点数量,a,开始点,b结束点 h=(b-a)/count_node;%h步长 np1=2*count_node+1;%包括端点,总共的2n+1点 x=a:h/2:b;%x自变量 c=ones(np1,1); for i=1:np1 if(mod(i,2
8、)=0)%偶数 c(i,1)=2/3; else%奇数 c(i,1)=1/3; end end c(1,1)=1/6;c(np1,1)=1/6;%两个端点都取1/6 syms symbol_x; fx=subs(f,symbol_x,x); result=vpa(h*fx*c,9);代码3，检验复合梯形公式和复合simpson公式function result1=test_integration() %测试复合梯形公式，复合simpson公式,用到trapezoid_integration,simpson_integration syms symbol_x; f1=sin(symbol_x)/
9、symbol_x; syms A1;%用于输出f1积分结果 for i=1:9 A1(i,1)=trapezoid_integration(f1,i,1e-9,1); A1(i,2)=simpson_integration(f1,i,1e-9,1); end k=10; for i=10:10:100 A1(k,1)=trapezoid_integration(f1,i,1e-9,1); A1(k,2)=simpson_integration(f1,i,1e-9,1); k=k+1; end IS(2)=int(f1,symbol_x,0,1); vpa(IS(2) result1=A1;结果
10、分析 从上面的计算结果可以看出复合梯形公式和复合simpson公式的稳定性，并且步长越小精度越高，当n趋于正无穷，即步长趋于0时，上述两公式的计算值都收敛到积分值。 从收敛的速度来看，复合simpson公式优于复合梯形公式。实验四 线方程组的直接解法一、问题提出：1、 三对角形线性方程组 二、要求 1、 对上述三个方程组分别利用Gauss顺序消去法与Gauss列主元消去法；平方根法与改进平方根法；追赶法求解（选择其一）； 2、 应用结构程序设计编出通用程序； 3、 比较计算结果，分析数值解误差的原因； 4、 尽可能利用相应模块输出系数矩阵的三角分解式。 三、目的和意义 1、通过该课题的实验，体
11、会模块化结构程序设计方法的优点； 2、运用所学的计算方法，解决各类线性方程组的直接算法； 3、提高分析和解决问题的能力，做到学以致用； 通过三对角形线性方程组的解法，体会稀疏线性方程组解法的特点。 四、实验步骤：高斯顺序消去法:function result=Gauss(A,d) %高斯顺序消元法 n,m=size(A) k=ones(1,n);%用于记录-A(j,i)/A(i,i) for i=1:n if(A(i,i)=0) 主元为0，不能使用高斯顺序消元法 return; end for j=i+1:n k(j)=-A(j,i)/A(i,i) A(j,:)=A(i,:).*k(j)+A(
12、j,:);%加到第j行 d(j)=d(i).*k(j)+d(j); end end x=zeros(n,1) x(n)=d(n)/A(n,n); for i=n-1:-1:1 x(i)=(d(i)-A(i,:)*x)/A(i,i); end result=xGauss列主元消去法function result=Gauss2(A,d) %高斯列主元消元法 n,m=size(A); d_length,d_width=size(d); if(d_length=n | d_width=1) 方程右端向量输入有误 return; end if(n=m) 系数矩阵不是方阵 return; end Aand
13、D=A d;%增广阵 k=ones(1,n);%用于记录-A(j,i)/A(i,i) record=1:n;%用于记录行的交换情况 for i=1:n %选取第i列的最大元进行换位 l,index=max(abs(AandD(i:n,i); if(index=i)temp_record=record(index); record(index)=record(i+index-1); record(i+index-1)=temp_record; temp=AandD(i,:); AandD(i,:)=AandD(i+index-1,:); AandD(i+index-1,:)=temp; end
14、%顺序消元 if(AandD(i,i)=0) 主元为0，不能使用高斯列主元顺序消元法 return; end for j=i+1:n k(j)=-AandD(j,i)/AandD(i,i); AandD(j,:)=AandD(i,:).*k(j)+AandD(j,:);%加到第j行 end endA=AandD(1:n,1:n); d=AandD(1:n,n+1); x=zeros(n,1); x(n)=d(n)/A(n,n); for i=n-1:-1:1 x(i)=(d(i)-A(i,:)*x)/A(i,i); end %由于不是顺序消元，所以还要再换回来 for i=1:n result
15、(i)=x(record(i); end追赶法function result=chase_after(A,d) n,m=size(A); d_length,d_width=size(d); if(d_length=n | d_width=1) 方程右端向量输入有误 return; end if(n=m) 系数矩阵不是方阵 return; end for i=2:ntoCheck1=diag(A,i); toCheck2=diag(A,-i); if(norm(toCheck1)=0|norm(toCheck2)=0) 系数矩阵不是三对角阵 return; end endb=diag(A); a
16、=zeros(n,1); a(2:n)=diag(A,-1); c=diag(A,1); if(b(1)=0) 因顺序主子式为0,不能进行Crout分解 return; endl(1)=b(1); y(1)=d(1)/l(1); u(1)=c(1)/l(1); for i=2:n l(i)=b(i)-a(i)*u(i-1); if(l(i)=0) 因顺序主子式为0,不能进行Crout分解 return; end y(i)=(d(i)-y(i-1)*a(i)/l(i); if(i=n) u(i)=c(i)/l(i); end end x(n)=y(n); for i=n-1:-1:1 x(i)=
17、y(i)-u(i)*x(i+1); end y=y;result=x;检验三种直接解线性方程组的方法function test_chase_after A= 4 -1 0 0 0 0 0 0 0 0; -1 4 -1 0 0 0 0 0 0 0; 0 -1 4 -1 0 0 0 0 0 0; 0 0 -1 4 -1 0 0 0 0 0; 0 0 0 -1 4 -1 0 0 0 0; 0 0 0 0 -1 4 -1 0 0 0; 0 0 0 0 0 -1 4 -1 0 0; 0 0 0 0 0 0 -1 4 -1 0; 0 0 0 0 0 0 0 -1 4 -1; 0 0 0 0 0 0 0 0
18、 -1 4; d=7 5 -13 2 6 -12 14 -4 5 -5;%按列存储 L,U = lu(A) %检验高斯顺序消去法 result_Gauss=Gauss(A,d) %检验高斯列主元消去法 result_Gauss2=Gauss2(A,d) %检验追赶法 result_chase=chase_after(A,d)运行结果：高斯顺序消去法的计算结果： result_Gauss = 2 1 -3 1.7849e-016 1 -2 3 -1.4874e-016 1 -1高斯列主元消去法的计算结果： result_Gauss = 2 1 -3 1.7849e-016 1 -2 3 -1.4
19、874e-016 1 -1追赶法的计算结果： result_Gauss = 2 1 -3 1.7849e-016 1 -2 3 -1.4874e-016 1 -1结果分析：高斯顺序消去法与精确解，比较其误差是1.7849*10(-16), 该误差可以忽略，事实上这是由计算机的除法引起的。而在这个问题中，因为顺序消元的除法中用到的数k的绝对值都小于1，这就避免使用绝对值较大的数，使得总体舍入误差得到有效控制，所以高斯顺序消元法的效果也很好。 同样高斯列主元消去法与精确解的误差都是（1）式的结果，说明在这个问题中误差都控制得很好。实验五 解线性方程组的迭代法一、问题提出 对实验四所列目的和意义的线
20、性方程组，试分别选用Jacobi 迭代法，Gauss-Seidel迭代法和SOR方法计算其解。 二、要求 1、体会迭代法求解线性方程组，并能与消去法做以比较； 2、分别对不同精度要求，如由迭代次数体会该迭代法的收敛快慢；3、对方程组2，3使用SOR方法时，选取松弛因子=0.8，0.9，1，1.1，1.2等，试看对算法收敛性的影响，并能找出你所选用的松弛因子的最佳者； 4、给出各种算法的设计程序和计算结果。 三、目的和意义 1、通过上机计算体会迭代法求解线性方程组的特点，并能和消去法比较； 2、运用所学的迭代法算法，解决各类线性方程组，编出算法程序； 3、体会上机计算时，终止步骤或k （给予的迭
21、代次数），对迭代法敛散性的意义； 4、 体会初始解，松弛因子的选取，对计算结果的影响。四、实验步骤：Jacobi迭代法function result=Jacobi(A,b,x0,yupsilong) %Jacobi迭代方法,x0是迭代用的初始向量,yupsilong是解的误差上限 n,m=size(A); b_length,b_width=size(b); if(b_length=n | b_width=1) 方程右端向量输入有误 return; end if(n=m) 系数矩阵不是方阵 return; end if(det(A)=0) 系数矩阵是奇异阵 return end %分解A=D-L
22、-U;D是对角元，L是对角元为0的下三角阵，U是对角元为0的上三角阵 L=tril(A,-1); U=triu(A,1); D=A-L-U; L=-L; U=-U; %计算迭代矩阵B0=D(-1); B=B0*(L+U); g=B0*b; x=x0; delta=1; count_iteration=0; while delta=yupsilong if(count_iteration1000) 超过迭代次数上限，认为算法此时不收敛 return; end x1=B*x+g; delta=norm(x1-x); x=x1; count_iteration=count_iteration+1;
23、end count_iteration result=x;Gauss_Seidol迭代法function result=Gauss_Seidol(A,b,x0,yupsilong) %Gauss_Seidol方法,x0是迭代用的初始向量,yupsilong是解的误差上限 n,m=size(A); b_length,b_width=size(b); if(b_length=n | b_width=1) 方程右端向量输入有误 return; end if(n=m) 系数矩阵不是方阵 return; end if(det(A)=0) 系数矩阵是奇异阵 return end %分解A=D-L-U;D是
24、对角元，L是对角元为0的下三角阵，U是对角元为0的上三角阵 L=tril(A,-1); U=triu(A,1); D=A-L-U; L=-L; U=-U;%计算迭代矩阵 B0=(D-L)(-1); B=B0*U; g=B0*b; x=x0; delta=1; count_iteration=0; while delta=yupsilong if(count_iteration1000) 超过迭代次数上限，认为算法此时不收敛 return; end x1=B*x+g; delta=norm(x1-x); x=x1; count_iteration=count_iteration+1; end c
25、ount_iteration result=x;SOR迭代法function result=SOR(A,b,x0,yupsilong,omiga) %SOR迭代方法,x0是迭代用的初始向量,yupsilong是解的误差上限,omiga是松弛因子 n,m=size(A); b_length,b_width=size(b); if(b_length=n | b_width=1) 方程右端向量输入有误 return; end if(n=m) 系数矩阵不是方阵 return; end if(det(A)=0) 系数矩阵是奇异阵 return end %分解A=D-L-U;D是对角元，L是对角元为0的下
26、三角阵，U是对角元为0的上三角阵 L=tril(A,-1); U=triu(A,1); D=A-L-U; L=-L;U=-U; %计算迭代矩阵 B0=(D-omiga*L)(-1); B=B0*(1-omiga)*D+omiga*U); g=omiga*B0*b; x=x0; delta=1; count_iteration=0; while delta=yupsilong if(count_iteration1000) 超过迭代次数上限，认为算法此时不收敛 return; end x1=B*x+g; delta=norm(x1-x); x=x1; count_iteration=count_
27、iteration+1; end count_iteration result=x;检验Jocobi,GS,SOR迭代方法function test_Gauss_Seidol clc; A= 4 -1 0 0 0 0 0 0 0 0; -1 4 -1 0 0 0 0 0 0 0; 0 -1 4 -1 0 0 0 0 0 0; 0 0 -1 4 -1 0 0 0 0 0; 0 0 0 -1 4 -1 0 0 0 0; 0 0 0 0 -1 4 -1 0 0 0; 0 0 0 0 0 -1 4 -1 0 0; 0 0 0 0 0 0 -1 4 -1 0; 0 0 0 0 0 0 0 -1 4 -1
28、; 0 0 0 0 0 0 0 0 -1 4; b=7 5 -13 2 6 -12 14 -4 5 -5; n,m=size(A); yupsilong=1e-3 1e-4 1e-5 1e-6 1e-7 1e-8 1e-9 1e-10; x_real=2 1 -3 0 1 -2 3 0 1 -1;%精确解 测试初始点对迭代的影响for i=1:8 x0=ones(n,1)*10(i-1); x0 result_Jacobi=(Jacobi(A,b,x0,yupsilong(3); if(norm(result_Jacobi-x_real)yupsilong(3) 误差过大 end result
29、_GS=(Gauss_Seidol(A,b,x0,yupsilong(3); if(norm(result_GS-x_real)yupsilong(3) 误差过大 end result_SOR=(SOR(A,b,x0,yupsilong(3),1.1); if(norm(result_GS-x_real)yupsilong(3) 误差过大 end end 测试yupsilong对迭代的影响x0=ones(n,1); for i=1:length(yupsilong) (yupsilong(i) result_Jacobi=(Jacobi(A,b,x0,yupsilong(i); result_GS=(Gauss_Seidol(A,b,x0,yupsilong
