数值分析报告资料

《数值分析》课程名称数值分析实验报告专业班级应用数学111班姓名陈伟学号8教学教师岳雪芝实验二函数逼近与曲线拟合一、问题提出从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量与时间t的拟合曲线。t(分)05101520253035404550554(10)y01.272.162.863.443.874.154.374.514.584.024.64二、要求1、用最小二乘法进行曲线拟合；2、近似解析表达式为23123()tatatat；3、打印出拟合函数()t，并打印出()jt与()jyt的误差，1,2,,12j；4、另外选取一个近似表达式，尝试拟合效果的比较；5、*绘制出曲线拟合图。三、目的和意义1、掌握曲线拟合的最小二乘法；2、最小二乘法亦可用于解超定线代数方程组；3、探索拟合函数的选择与拟合精度间的关系。四、实验步骤:输入代码：t=[0510152025303540455055];y=[01.272.162.863.443.874.154.374.514.584.024.64];fori=1:12A(i,1)=t(i);A(i,2)=t(i).^2;A(i,3)=t(i).^3;endp_star=A\y';y_star=p_star(1)*t+p_star(2)*t.^2+p_star(3)*t.^3;s_star=0;fori=1:12s_star=s_star+(y_star(i)-y(i))^2;endplot(t,y,'*',t,y_star,'g');比较拟和结果，分别取二次，三次和四次多项式，再做最小二乘法，代码如下：t=[0510152025303540455055];y=[01.272.162.863.443.874.154.374.514.584.024.64];%3次拟合p3=polyfit(t,y,3)y_new=polyval(p3,t);s_3=0;fori=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;fori=1:12s_2=s_2+(y_new2(i)-y(i))^2;end%4次拟合p4=polyfit(t,y,4);y_new4=polyval(p4,t);s_4=0;fori=1:12s_4=s_4+(y_new4(i)-y(i))^2;end比较4种拟合函数，结论：常数项为0的三次拟合函数在保证拟合精度的同时，保证0点的函数值仍为0，是四种拟合方法中最好的一种。原图与曲线拟合图如下所示：0102030405060-10123456x10-4分析:从上面的拟合中也可以知道多项式拟合误差平方和随着拟合多项式次数的增加而逐渐减少，拟合曲线更靠近实际数据，拟合更准确。实验三数值积分与数值微分一、基本题选用复合梯形公式，复合Simpson公式，Romberg算法，计算（1）10sinI=((0)1,0.9460831)xdxfIx二、要求1、编制数值积分算法的程序；2、分别用两种算法计算同一个积分，并比较其结果；3、分别取不同步长()/hban，试比较计算结果（如n=10,20等）；4、给定精度要求ε，试用变步长算法，确定最佳步长。三、目的和意义1、深刻认识数值积分法的意义；2、明确数值积分精度与步长的关系；3、根据定积分的计算方法，结合专业考虑给出一个二重积分的计算问题。四、实验步骤：代码1，复合梯形公式functionresult=trapezoid_integration(f,count_node,a,b)%compositetrapezoidintegration,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;%两个端点都取一半symssymbol_x;fx=subs(f,symbol_x,x);result=vpa(h*fx*c,9);代码2，复合simpson公式functionresult=simpson_integration(f,count_node,a,b)%compositesimpsonintegration,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);fori=1:np1if(mod(i,2)==0)%偶数c(i,1)=2/3;else%奇数c(i,1)=1/3;endendc(1,1)=1/6;c(np1,1)=1/6;%两个端点都取1/6symssymbol_x;fx=subs(f,symbol_x,x);result=vpa(h*fx*c,9);代码3，检验复合梯形公式和复合simpson公式function[result1]=test_integration()%测试复合梯形公式，复合simpson公式,用到trapezoid_integration,simpson_integrationsymssymbol_x;f1=sin(symbol_x)/symbol_x;symsA1;%用于输出f1积分结果fori=1:9A1(i,1)=trapezoid_integration(f1,i,1e-9,1);A1(i,2)=simpson_integration(f1,i,1e-9,1);endk=10;fori=10:10:100A1(k,1)=trapezoid_integration(f1,i,1e-9,1);A1(k,2)=simpson_integration(f1,i,1e-9,1);k=k+1;endIS(2)=int(f1,symbol_x,0,1);vpa(IS(2))result1=A1;结果分析从上面的计算结果可以看出复合梯形公式和复合simpson公式的稳定性，并且步长越小精度越高，当n趋于正无穷，即步长趋于0时，上述两公式的计算值都收敛到积分值。从收敛的速度来看，复合simpson公式优于复合梯形公式。实验四线方程组的直接解法一、问题提出：1、三对角形线性方程组123456789104100000000141000000001410000000014100000000141000000001410000000014100000000141000000001410000000014xxxxxxxxxx7513261214455*(2,1,3,0,1,2,3,0,1,1)Tx二、要求1、对上述三个方程组分别利用Gauss顺序消去法与Gauss列主元消去法；平方根法与改进平方根法；追赶法求解（选择其一）；2、应用结构程序设计编出通用程序；3、比较计算结果，分析数值解误差的原因；4、尽可能利用相应模块输出系数矩阵的三角分解式。三、目的和意义1、通过该课题的实验，体会模块化结构程序设计方法的优点；2、运用所学的计算方法，解决各类线性方程组的直接算法；3、提高分析和解决问题的能力，做到学以致用；通过三对角形线性方程组的解法，体会稀疏线性方程组解法的特点。四、实验步骤：高斯顺序消去法:functionresult=Gauss(A,d)%高斯顺序消元法[n,m]=size(A)k=ones(1,n);%用于记录-A(j,i)/A(i,i)fori=1:nif(A(i,i)==0)'主元为0，不能使用高斯顺序消元法'return;endforj=i+1:nk(j)=-A(j,i)/A(i,i)A(j,:)=A(i,:).*k(j)+A(j,:);%加到第j行d(j)=d(i).*k(j)+d(j);endendx=zeros(n,1)x(n)=d(n)/A(n,n);fori=n-1:-1:1x(i)=(d(i)-A(i,:)*x)/A(i,i);endresult=xGauss列主元消去法functionresult=Gauss2(A,d)%高斯列主元消元法[n,m]=size(A);[d_length,d_width]=size(d);if(d_length~=n||d_width~=1)'方程右端向量输入有误'return;endif(n~=m)'系数矩阵不是方阵'return;endAandD=[Ad];%增广阵k=ones(1,n);%用于记录-A(j,i)/A(i,i)record=1:n;%用于记录行的交换情况fori=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%顺序消元if(AandD(i,i)==0)'主元为0，不能使用高斯列主元顺序消元法'return;endforj=i+1:nk(j)=-AandD(j,i)/AandD(i,i);AandD(j,:)=AandD(i,:).*k(j)+AandD(j,:);%加到第j行endendA=AandD(1:n,1:n);d=AandD(1:n,n+1);x=zeros(n,1);x(n)=d(n)/A(n,n);fori=n-1:-1:1x(i)=(d(i)-A(i,:)*x)/A(i,i);end%由于不是顺序消元，所以还要再换回来fori=1:nresult(i)=x(record(i));end追赶法functionresult=chase_after(A,d)[n,m]=size(A);[d_length,d_width]=size(d);if(d_length~=n||d_width~=1)'方程右端向量输入有误'return;endif(n~=m)'系数矩阵不是方阵'return;endfori=2:ntoCheck1=diag(A,i);toCheck2=diag(A,-i);if(norm(toCheck1)~=0||norm(toCheck2)~=0)'系数矩阵不是三对角阵'return;endendb=diag(A);a=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);fori=2:nl(i)=b(i)-a(i)*u(i-1);if(l(i)==0)'因顺序主子式为0,不能进行Crout分解'return;endy(i)=(d(i)-y(i-1)*a(i))/l(i);if(i~=n)u(i)=c(i)/l(i);endendx(n)=y(n);fori=n-1:-1:1x(i)=y(i)-u(i)*x(i+1);endy=y';result=x';检验三种直接解线性方程组的方法functiontest_chase_afterA=