数值稳定性验证实验报告

实验课程：数值计算方法专业：数学与应用数学班级：08070141学号：37*名：***中北大学理学院实验1 赛德尔迭代法【实验目的】熟悉用塞德尔迭代法解线性方程组 【实验内容】1.了解MATLAB 语言的用法2.用塞德尔迭代法解下列线性方程组1234123412341234541012581034x x x x x x x x x x x x x x x x ---=-⎧⎪-+--=⎪⎨--+-=⎪⎪---+=⎩ 【实验所使用的仪器设备与软件平台】计算机，MATLAB7.0【实验方法与步骤】1.先找出系数矩阵A ，将前面没有算过的x j 分别和矩阵的(,)A i j 相乘,然后将累加的和赋值给sum ，即(),j sum sum A i j x =+⋅.算出()/(,)i i x b sum A i i =-,依次循环,算出所有的i x 。

2.若i x 前后两次之差的绝对值小于所给的误差限ε，则输出i x .否则重复以上过程，直到满足误差条件为止.【实验结果】(A 是系数矩阵，b 是右边向量，x 是迭代初值，ep 是误差限)function y=seidel(A,b,x,ep) n=length(b); er=1; k=0; while er>=epk=k+1;for i=[1:1:n]q=x(i);sum=0;for j=[1:1:n]if j~=isum=sum+A(i,j)*x(j);endendx(i)=(b(i)-sum)/A(i,i);er=abs(q-x(i));endendfprintf('迭代次数k=%d\n',k)disp(x')【结果分析与讨论】>> A=[5 -1 -1 -1;-1 10 -1 -1;-1 -1 5 -1;-1 -1 -1 10]; b=[-4 12 8 34];seidel(A,b,[0 0 0 0],1e-3)迭代次数k=60.998978494300021.999584568676492.999531397434353.99980944604109实验课程：数值计算方法专业：数学与应用数学班级：08070141学号：37姓名：汪鹏飞中北大学理学院实验2 最小二乘法的拟合【实验目的】熟悉最小二乘法的拟合方法 【实验内容】1.了解MATLAB 语言的用法2.给出数据希望用一次，二次多项式利用最小二乘法拟合这些数据，试写出正规方程组，并求出最小平方逼近多项式。

【实验所使用的仪器设备与软件平台】计算机，MATLAB7.0，【实验方法与步骤】1.先算出 1n j j i i l x ==∑，再算出 (1)1nj j i i i b x y -==∑2.在正规方程组20121111231012111111201211111n n n nmi i m i ii i i i n n n n nm i i i m i i ii i i i i n n n n nm m m m m m i i i m i i ii i i i i a n a x a x a x y a x a x a x a x x y a x a x a x a x x y ====+=====+++=====⎧++++=⎪⎪⎪++++=⎪⎨⎪⎪⎪++++=⎪⎩∑∑∑∑∑∑∑∑∑∑∑∑∑∑ 中令1n j j i i l x ==∑，(1)1nj j i i i b x y -==∑3.令正规方程组的系数矩阵为A ，当j m ≤时，将j l 的值赋给(1,1)A j +， 当j>m 时，将j l 的值赋给(1,1)A j m m +-+，再将矩阵A 的其他元素写出来，于是正规矩阵可写成Aa b =,最后用1a A b -=即可算出向量a,向量a 的元素依次是常数项,一次项的系数,二次项的系数m次项的系数.4.由系数即可写出拟合的多项式曲线.【实验结果】（x，y为给定的对应值，m是要求的拟合次数)function a=zxecf(x,y,m)n=length(x);A(1,1)=n;for j=[1:1:2*m]l(j)=0;for i=[1:1:n]l(j)=l(j)+x(i)^j;endif j<=mA(1,1+j)=l(j);else A(j+1-m,m+1)=l(j);endendfor j=[1:1:m+1]b(j)=0;for i=[1:1:n]b(j)=b(j)+y(i)*x(i)^(j-1);endendfor i=[2:1:m+1]k=-1;for j=[m+1:-1:1]k=k+1;A(i,j)=l(m+i-1-k);end end a=A\b';【结果分析与讨论】>> x=[-1.00 -0.75 -0.50 -0.25 0 0.25 0.50 0.75 1.00];y=[-0.2209 0.3295 0.8826 1.4392 2.0003 2.5645 3.1334 3.7061 4.2836]; a1=zxecf(x,y,1) a1 =2.01314444444444 2.25164666666667>> a2=zxecf(x,y,2) a2 =2.00010043290043 2.25164666666667 0.03130562770563拟合的一次多项式是() 2.0131 2.2516y x x =+拟合的二次多项式是2() 2.0001 2.25160.0313y x x x =++优缺点:该方法可以通过改变m 的值来将数据拟合成任意次数的多项式函数。

实验课程：数值计算方法专业：数学与应用数学班级：08070141学号：37姓名：汪鹏飞中北大学理学院实验3 龙贝格算法求积分【实验目的】熟悉龙贝格方法求积分 【实验内容】1.了解MATLAB 语言的用法2.用龙贝格方法求积分21-5(10)x I edx -=⎰要求误差不超过【实验所使用的仪器设备与软件平台】计算机，MATLAB7.0 【实验方法与步骤】1.计算()f a 和()f b ，算出1T .2.将区间[],a b 分半，计算2a b f +⎛⎫⎪⎝⎭，算出2T 和1S .3.再将子区间分半，计算 4b a f a -⎛⎫+ ⎪⎝⎭，34b a fa -⎛⎫+⋅⎪⎝⎭，算出4T ,2S 和1C . 4.再将子区间分半，计算 8b a f a -⎛⎫+ ⎪⎝⎭,38b a f a -⎛⎫+⋅ ⎪⎝⎭,58b a f a -⎛⎫+⋅ ⎪⎝⎭,78b a f a -⎛⎫+⋅ ⎪⎝⎭，算出8T , 4S , 2C 和1R .5.再将子区间分半,算出16T , 8S , 4C 和 2R ,不断将子区间分半，重复以 上过程，算出4R , 8R , 16R ，一直到相邻两个R 值之差的绝对值不超过给定的误差ε为止，最后一次算出的R 值即为所求。

【实验结果】(a 是下限，b 是上限，eps 是误差限)定义函数function y=intf(x)y=2*(exp(-x^2))/(pi^0.5);主程序function J=romberg(a,b,eps)er=1;k=3;h=b-a;T(1)=h*(intf(a)+intf(b))/2;h1=h/2;T(2)=T(1)/2+h1*intf(a+h1);S(1)=4/(4-1)*T(2)-1/(4-1)*T(1);h2=h1/2;T(4)=T(2)/2+h2*(intf(a+h2)+intf(a+3*h2));S(2)=4/(4-1)*T(4)-1/(4-1)*T(2);C(1)=4^2/(4^2-1)*S(2)-1/(4^2-1)*S(1);h3=h2/2;T(8)=T(4)/2+h3*(intf(a+h3)+intf(a+3*h3)+intf(a+5*h3)+intf(a+7*h3)); S(4)=4/(4-1)*T(8)-1/(4-1)*T(4);C(2)=4^2/(4^2-1)*S(4)-1/(4^2-1)*S(2);R(1)=4^3/(4^3-1)*C(2)-1/(4^3-1)*C(1);while er>=epsH=0;k=k+1;u=h/2^k;x=a+u;i=0;while x<bH=H+intf(x);i=i+1;x=a+u*(2*i+1);endT(2^k)=T(2^(k-1))/2+h*H/2^k;S(2^(k-1))=4/(4-1)*T(2^k)-1/(4-1)*T(2^(k-1));C(2^(k-2))=4^2/(4^2-1)*S(2^(k-1))-1/(4^2-1)*S(2^(k-2));R(2^(k-3))=4^3/(4^3-1)*C(2^(k-2))-1/(4^3-1)*C(2^(k-3));er=abs(R(2^(k-3))-R(2^(k-4)));endI=R(2^(k-3));fprintf('I=%10.8f\n',I)【结果分析与讨论】>> romberg(0,1,1e-4)I=0.84270079实验课程：数值计算方法专业：数学与应用数学班级：08070141学号：37姓名：汪鹏飞中北大学理学院实验4 标准四阶R-K 方法【实验目的】熟悉四阶R-K 方法【实验内容】1.了解MATLAB 语言的用法2.取步长0.2h =，用标准四阶R-K 方法解初值问题3,011(0)1y y x x y ⎧'=≤≤⎪+⎨⎪=⎩ 【实验所使用的仪器设备与软件平台】计算机，MATLAB7.0【实验方法与步骤】1.先将积分区除以步长确定分点个数b a n h-=. 2.利用标准四阶R-K 公式，算出()1,n n k h f x y =⋅，12,22n n k h k h f x y ⎛⎫=⋅++ ⎪⎝⎭; 23,22n n k h k h f x y ⎛⎫=⋅++ ⎪⎝⎭， ()43,n n k h f x h y k =⋅++.3.利用()112341226n n y y k k k k +=+⋅+++算出下一个节点处的y 值. 4.重复以上步骤，算出每个分点处y 的值.5.输出每个分点处y 的值.【实验结果】(a是下限，b是上限，h是步长，y0是初值) 定义函数function f=func(x,y)f=(3*y)/(1+x);主程序function q=LK(a,b,h,y0)n=(b-a)/h;fprintf('a=%6.4f,y0=%10.8f\n',a,y0)for i=[0:1:n-1]k1=h*func(a,y0);k2=h*func(a+h/2,y0+k1/2);k3=h*func(a+h/2,y0+k2/2);k4=h*func(a+h,y0+k3);y0=y0+1/6*(k1+2*k2+2*k3+k4);fprintf('a=%6.4f,y0=%10.8f\n',a+h,y0)a=a+h;end【结果分析与讨论】>> LK(0,1,0.2,1)a=0.0000,y0=1.00000000a=0.2000,y0=1.72754821a=0.4000,y0=2.74295130a=0.6000,y0=4.09418136a=0.8000,y0=5.82921073a=1.0000,y0=7.99601214。