数值计算方法实验报告

【实验目的】

熟悉用塞德尔迭代法解线性方程组 【实验内容】

1.了解MATLAB语言的用法

2.用塞德尔迭代法解下列线性方程组

5x1x2x3x4x10xxx1234x1x25x3x4x1x2x310x4412 834【实验所使用的仪器设备与软件平台】

计算机，MATLAB7.0

【实验方法与步骤】

1.先找出系数矩阵A，将前面没有算过的xj分别和矩阵的A(i,j)相乘,然后将累加的和赋值给sum，即sumsumAi,jxj.算出

xi(bisum)/A(i,i),依次循环,算出所有的xi 。

2.若xi前后两次之差的绝对值小于所给的误差限，则输出xi.否则重复以上过程，直到满足误差条件为止.

【实验结果】

(A是系数矩阵，b是右边向量，x是迭代初值，ep是误差限)

function y=seidel(A,b,x,ep) n=length(b); er=1; k=0; while er>=ep

k=k+1; for i=[1:1:n] q=x(i); sum=0; for j=[1:1:n] if j~=i

sum=sum+A(i,j)*x(j); end end

x(i)=(b(i)-sum)/A(i,i); er=abs(q-x(i)); end end

fprintf('迭代次数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=6 0.99897849430002 1.99958456867649 2.99953139743435 3.99980944604109

实验课程： 数值计算方法 专 业： 数学与应用数学 班 级： 08070141 学 号： 37 姓 名： 汪鹏飞

中北大学理学院

实验2 最小二乘法的拟合

【实验目的】

熟悉最小二乘法的拟合方法 【实验内容】

1.了解MATLAB语言的用法 2.给出数据 x y -1.00 -0.2209 -0.75 -0.50 -0.25 0 0.25 0.50 0.75 1.00 0.3295 0.8826 1.4392 2.0003 2.5645 3.1334 3.7061 4.2836 希望用一次，二次多项式利用最小二乘法拟合这些数据，试写出正规方程组，并求出最小平方逼近多项式。

【实验所使用的仪器设备与软件平台】

计算机，MATLAB7.0，

【实验方法与步骤】

1.先算出 ljxi，再算出 bjxi(j1)yi

ji1i1nn 2.在正规方程组

nnnn2ma0na1xia2xiamxiyii1i1i1i1nnnnn23m1a0xia1xia2xiamxixiyi  i1i1i1i1i1nnnnnmm1m2mmma0xia1xia2xiamxixiyii1i1i1i1i1 中令ljxi，bjxi(j1)yi

ji1i1nn 3.令正规方程组的系数矩阵为A，当jm时，将lj的值赋给A(1,1j)， 当j>m时，将lj的值赋给A(j1m,m1)，再将矩阵A的其他元素写出来，于是正规矩阵可写成Aab,最后用aA1b即可算出向量a,向量a的

元素依次是常数项,一次项的系数,二次项的系数4.由系数即可写出拟合的多项式曲线.

m次项的系数.

【实验结果】

（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; end if j<=m

A(1,1+j)=l(j); else A(j+1-m,m+1)=l(j); end end

for j=[1:1:m+1] b(j)=0; for i=[1:1:n]

b(j)=b(j)+y(i)*x(i)^(j-1); end end

for 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

拟合的一次多项式是y(x)2.01312.2516x

拟合的二次多项式是y(x)2.00012.2516x0.0313x2

优缺点:

该方法可以通过改变m的值来将数据拟合成任意次数的多项式函数。

实验课程： 数值计算方法 专 业： 数学与应用数学 班 级： 08070141 学 号： 37 姓 名： 汪鹏飞

中北大学理学院

实验3 龙贝格算法求积分

【实验目的】

熟悉龙贝格方法求积分 【实验内容】

1.了解MATLAB语言的用法 2.用龙贝格方法求积分

I21x20edx(要求误差不超过10-5)

【实验所使用的仪器设备与软件平台】

计算机，MATLAB7.0 【实验方法与步骤】

1.计算fa和fb，算出T1.

ab 2.将区间a,b分半，计算f，算出T2和S1.

2ba 3.再将子区间分半，计算 fa，f4ba 4.再将子区间分半，计算 fa,

8baa3，算出T4,S2和C1. 4bafa5,

8bafa3,

8bafa7，算出T8, S4, C2 和R1.

8 5.再将子区间分半,算出T16, S8, C4 和 R2,不断将子区间分半，重复以 上过程，算出R4, R8, R16，一直到相邻两个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>=eps H=0; k=k+1; u=h/2^k; x=a+u; i=0; while xH=H+intf(x); i=i+1; x=a+u*(2*i+1); end

T(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))); end

I=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.取步长h0.2，用标准四阶R-K方法解初值问题

3y,0x1y 1xy(0)1【实验所使用的仪器设备与软件平台】 计算机，MATLAB7.0 【实验方法与步骤】

1.先将积分区除以步长确定分点个数nba. h 2.利用标准四阶R-K公式，算出

k1hfxn,yn，

khk2hfxn,yn1;

22khk3hfxn,yn2，

22k4hfxnh,ynk3.

1 3.利用yn1ynk12k22k3k4算出下一个节点处的y值.

6 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.00000000 a=0.2000,y0=1.72754821 a=0.4000,y0=2.74295130 a=0.6000,y0=4.09418136 a=0.8000,y0=5.82921073 a=1.0000,y0=7.99601214