西南交通大学
数值分析2015上机实习报告
2015年11月
数值分析2015上机实习报告
目录
第2题 (1)
1. 程序 (1)
2. 结果分析 (2)
第3题 (5)
1. 程序 (6)
2. 结果分析 (7)
第5题 (8)
1. 程序 (8)
2. 结果分析 (9)
第2题
2. 某过程测涉及两变量x 和y, 拟分别用插值多项式和多项式拟合给出其对应规律的近似多项式，已知xi与yi之间的对应数据如下，xi=1,2,…,10
yi = 34.6588 40.3719 14.6448 -14.2721 -13.3570 24.8234 75.2795
103.5743 97.4847 78.2392
（1）请用次数分别为3，4，5，6的多项式拟合并给出最好近似结果f(x)。
（2）请用插值多项式给出最好近似结果
下列数据为另外的对照记录，它们可以作为近似函数的评价参考数据。
xi =
Columns 1 through 7
1.5000 1.9000
2.3000 2.7000
3.1000 3.5000 3.9000
Columns 8 through 14
4.3000 4.7000
5.1000 5.5000 5.9000
6.3000 6.7000
Columns 15 through 17
7.1000 7.5000 7.9000
yi =
Columns 1 through 7
42.1498 41.4620 35.1182 24.3852 11.2732 -1.7813 -12.3006
Columns 8 through 14
-18.1566 -17.9069 -11.0226 2.0284 19.8549 40.3626 61.0840
Columns 15 through 17
79.5688 93.7700 102.3677
1. 程序
（1）多项式拟合程序
n=input('输入所要拟合的阶数n=');
hold on;
x=1:10;
y=[34.6588 40.3719 14.6448 -14.2721 -13.3570 24.8234 75.2795 103.5743 97.4847 78.2392];
P=polyfit(x,y,n)
xi=1:.2:10;
yi=polyval(P,xi);
plot(xi,yi,x,y,'r*');
legend('拟合曲线','原始数据')
（2）拉格朗日插值多项式拟合程序
clc;
x=1:10;
y=[34.6588 40.3719 14.6448 -14.2721 -13.3570 24.8234 75.2795 103.5743 97.4847 78.2392]; plot(x,y,'r*');
hold on;
syms t;
n=length(x);
f = 0.0;
for i =1:n
l = y(i);
for j = 1:i-1
l = l*(t-x(j))/(x(i)-x(j));
end;
for j = i+1:n
l = l*(t-x(j))/(x(i)-x(j));
end;
f = f + l;
simplify(f);
f = collect(f);
f = vpa(f,6);
end
ti=1.0:0.4:10;
f=subs(f,'t',ti);
plot(ti,f)
legend('拟合曲线','原始数据')
title '拉格朗日插值多项式拟合'
2. 结果分析
（1）请用次数分别为3，4，5，6的多项式拟合并给出最好近似结果f(x)。
图2-1 3次多项式拟合结果
图2-2 4次多项式拟合结果
图2-3 5次多项式拟合结果
图2-4 6次多项式拟合结果
从绘制的图形来看，当采用6次多项式拟合的时候，拟合的曲线已经与所给出的点非常逼近了。
6次多项式拟合曲线为：f(x)=0.0194x6-0.5408x5+5.1137x4-16.8973x3-0.8670x2+66.3750x-18.6991（2）拉格朗日插值多项式给出最好近似结果
图1-4 拉格朗日插值多项式拟合
第3题
3．用雅格比法与高斯－赛德尔迭代法解下列方程组Ax=b，研究其收敛性，上机验证理论分析是否正确，比较它们的收敛速度，观察右端项对迭代收敛有无影响。
(1)A行分别为A1=[6,2,-1],A2=[1,4,-2],A3=[-3,1,4]；b1=[-3,2,4]T, b2=[100,-200,345]T,
(2) A行分别为A1=[1,0,8,0.8],A2=[0.8,1,0.8],A3=[0.8,0.8,1]；b1=[3,2,1]T, b2=[5,0,-10]T,
(3)A行分别为A1=[1,3],A2=[-7,1]；b=[4,6]T, 1. 程序
（1）雅格比法程序
function jacobi( )
clc;
clear;
A=[6 2 -1;1 4 -2;-3 1 4];
B=[-3;2;4];
Err_user=0.01;
N=500;
[m,n]=size(A);
X=zeros(n,1);
k=1;
while k<=N;
Xk=X;
for i=1:n
for j=1:n
if i~=j
AX(j)=A(i,j)*Xk(j);
end
end
Sum_AX=sum(AX);
AX=0;
X(i)=(B(i)-Sum_AX)/A(i,i);
end
E=max(abs(Xk-X));
if E break; end k=k+1; end disp(X); %显示迭代结果 disp(k); %显示迭代次数 （2）高斯－赛德尔迭代法程序 function GS( ) clc; clear; A=[6 2 -1;1 4 -2;-3 1 4]; B=[-3;2;4]; A=[1 0.8 0.8;0.8 1 0.8;0.8 0.8 1]; B1=[3;2;1]; B2=[5;0;-10]; Err_user=0.0001; N=50; [m,n]=size(A); X=zeros(n,1); k=1; while k<=N Xk=X; for i=1:n for j=1:n if i~=j AX(j)=A(i,j)*X(j); %ó? Jacobi ·¨?÷òa??±e end end sum_AX=sum(AX); AX=0;X(i)=(B1(i)-sum_AX)/A(i,i); end Er=max(abs(Xk-X)); if Er break; end k=k+1; end disp(X); disp(k); end 2. 结果分析 （1） 1）A行分别为A1=[6,2,-1],A2=[1,4,-2],A3=[-3,1,4]；b1=[-3,2,4]T 雅克比迭代：计算结果x=[-0.6363 0.5960 0.3737]T,.迭代次数16次。 高斯赛德迭代：计算结果x=[-0.6363 0.5959 0.3738]T,.迭代次数10次。 2）A行分别为A1=[1,0,8,0.8],A2=[0.8,1,0.8],A3=[0.8,0.8,1]；b1=[3,2,1]T 雅克比迭代：计算不收敛 高斯赛德迭代：计算结果x=[5.7690 0.7694 -4.2307]T,.迭代次数30次。 （2） A行分别为A1=[1,3],A2=[-7,1]；b1=[4,6]T。 雅克比迭代：计算结果不收敛 高斯赛德迭代：计算结果不收敛 通过对雅克比法和高斯法的上机编程实习，分析对比实验结果可得：在方程组Ax=b中，右端项对迭代收敛是有影响的，即当b增大时，迭代次数增加，收敛速度降低。并且通过对比可知，在相同条件下，高斯－赛德尔迭代法比雅克比迭代法收敛速度快；方法的选择也很重要，比如第二问，用雅克比迭代法是发散的，而用高斯迭代法则是收敛的。通过上机验证，我们也得出结论：理论分析是正确的，即当迭代矩阵的谱半径小于1时，迭代法是收敛的，而当迭代矩阵谱半径大于1的时，迭代法都是发散的。 第5题 5.用Ru n ge-Kutt a4阶算法对初值问题y/=-20*y，y(0)=1按不同步长求解，用于观察稳定区 间的作用，推荐两种步长h=0.1,0.2。 注：此方程的精确解为：y=e-20x 1.程序 %%%%%%%%%%%%%%%%%%%%% %Runge-Kutta4阶算法 %f=-20y %y（0）=1 %%%%%%%%%%%%%%%%%%%%% clc; clear; N=10;%设定节点个数 h=0.05;%设定步长 x(1)=0;%x0=0 y(1)=1;%y（0）=1 yy(1)=exp(-20*x(1));%y（0）的精确解 %%%%%%%%%%%%%%%%%%%%%%%开始用runge-kutta法计算 for i = 2:N K1 = -20*(y(i-1));%(xi,yi)点的导数为f=-20*y K2 = -20*(y(i-1)+K1*h/2); K3 = -20*(y(i-1)+K2*h/2); K4 = -20*(y(i-1)+K3*h); delta = h*(K1 + 2*K2 + 2*K3 + K4)/6; y(i) = y(i-1) + delta;%计算y（i）值 x(i) = x(i-1) + h;%计算下个节点的x（i)值 yy(i) = exp(-20*x(i));%计算y(i)的精确值 end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%存储计算结果 dlmwrite('f:\计算结果.xls', x, 'delimiter', '\t', 'precision', 8); dlmwrite('f:\计算结果.xls', y,'-append', 'delimiter', '\t', 'precision', 8); dlmwrite('f:\计算结果.xls', yy,'-append', 'delimiter', '\t', 'precision', 8); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%画图功能plot(x,y,'o',x,yy,'*'); xlabel('x'); ylabel('y'); legend('四阶Runge-Kutta法','精确解'); pause; close; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2. 结果分析 （1）步长h = 0.1 时结果 当步长取为0.1时，，计算结果如图1.2-2： 图5-1 h=0.05时计算结果 计算结果整理如表5-1所示： 表5-1 h=0.1时结果 ○21.2.3步长h = 0.2 时结果 当步长取为0.2时，，计算结果如图5-2： 图5-2 h=0.1时计算结果 计算结果整理如表5-2所示： 图5-2 h=0.2时计算结果 通过以上对Runge-Kutta法的应用，计算结果表明步长h的取值会影响算法的收敛性和稳定性： （1）Runge-Kutta法的步长h越长，计算结果的精确度越低，甚至计算结果不收敛； （2）当步长h在稳定区间时，误差逐步衰减。
相关主题
文本预览