数值分析上机实验报告总结归纳.docx
《数值分析上机实验报告总结归纳.docx》由会员分享,可在线阅读,更多相关《数值分析上机实验报告总结归纳.docx(25页珍藏版)》请在冰豆网上搜索。
数值分析上机实验报告总结归纳
实验报告一
题目:
非线性方程求解
摘要:
非线性方程的解析解通常很难给出,因此线性方程的数值解法就尤为重要。
本实验采用两种常见的求解方法二分法和Newton法及改进的Newton法。
前言:
(目的和意义)
掌握二分法与Newton法的基本原理和应用。
数学原理:
对于一个非线性方程的数值解法很多。
在此介绍两种最常见的方法:
二分法和Newton法。
对于二分法,其数学实质就是说对于给定的待求解的方程f(x),其在[a,b]上连续,f(a)f(b)<0,且f(x)在[a,b]内仅有一个实根x*,取区间中点c,若,则c恰为其根,否则根据f(a)f(c)<0是否成立判断根在区间[a,c]和[c,b]中的哪一个,从而得出新区间,仍称为[a,b]。
重复运行计算,直至满足精度为止。
这就是二分法的计算思想。
Newton法通常预先要给出一个猜测初值xo,然后根据其迭代公式产生逼近解x*的迭代数列{xk},这就是Newton法的思想。
当xo接近x*时收敛很快,但是当xo选择不好时,可能会发散,因此初值的选取很重要。
另外,若将该迭代公式改进为其中r为要求的方程的根的重数,这就是改进的Newton法,当求解已知重数的方程的根时,
在同种条件下其收敛速度要比Newton法快的多。
程序设计:
本实验采用Matlab的M文件编写。
其中待求解的方程写成function的方式,如下
functiony=f(x);
y=-x*x-sin(x);
写成如上形式即可,下面给出主程序。
二分法源程序:
clear
%%%给定求解区间
b=;
a=0;
%%%误差
R=1;
k=0;%迭代次数初值
while(R>5e-6);
c=(a+b)/2;
iff12(a)*f12(c)>0;
a=c;
else
b=c;
end
R=b-a;%求出误差
k=k+1;
end
x=c%给出解
Newton法及改进的Newton法源程序:
clear
%%%%输入函数
f=input('请输入需要求解函数>>','s')
%%%求解f(x)的导数df=diff(f);
%%%改进常数或重根数
miu=2;
%%%初始值x0
xO=input('inputinitialvaluexO>>');
k=0;%迭代次数
max=100;%最大迭代次数
R=eval(subs(f,'xO','x'));%求解f(xO),以确定初值xO时否就是解while(abs(R)>1e-8)
x1=xO-miu*eval(subs(f,'xO','x'))/eval(subs(df,'xO','x'));
R=x1-x0;
x0=x1;
k=k+1;
if(eval(subs(f,'xO','x'))<1e-1O);
break
end
ifk>max;%如果迭代次数大于给定值,认为迭代不收敛,重新输入初值ss=input('mayberesultiserror,chooseanewxO,y/n?
>>','s');
ifstrcmp(ss,'y')
xO=input('inputinitialvaluex0>>');
k=0;
else
break
end
end
end
k;%给出迭代次数x=xO;%给出解结果分析和讨论:
X2
1.用二分法计算方程sinx20在[1,2]内的根。
(5*106,下同)
计算结果为
x=;
f(x)=;
k=18;
由f(x)知结果满足要求,但迭代次数比较多,方法收敛速度比较慢。
2.用二分法计算方程x3x10在[1,]内的根。
计算结果为
x=;
f(x)=;
k=17;
由f(x)知结果满足要求,但迭代次数还是比较多。
3.用Newton法求解下列方程
a)xex10xo=;
计算结果为
x=;
f(x)=;k=4;
由f(x)知结果满足要求,而且又迭代次数只有4次看出收敛速度很快。
3
b)xx10xo=1;
2
c)(x1)(2x1)0xo=,xo=;
当灭0=时,计算结果为
x=;
f(x)=;k=4;
由f(x)知结果满足要求,而且又迭代次数只有4次看出收敛速度很快,实际上该方程确实有真解x=。
当灭0=时,计算结果为
x=;
f(x)=0;k=9;
由f(x)知结果满足要求,实际上该方程确实有真解x=,但迭代次数增多,实际上当取X0〉时,x~1就变成了方程的另一个解,这说明Newton法收敛与初值很有关系,有的时候甚至可能
不收敛。
4.用改进的Newton法求解,有2重根,取2
2
(x1)(2x1)0X0=;并与3.中的c)比较结果。
当X0=时,程序死循环,无法计算,也就是说不收敛。
改1.5时,结果收敛为
x=;
f(x)=;k=16;
显然这个结果不是很好,而且也不是收敛至方程的2重根上。
当X0=时,结果收敛为
x=;
f(x)=;k=4;
这次达到了预期的结果,这说明初值的选取很重要,直接关系到方法的收敛性,实际上直接用Newton法,在给定同样的条件和精度要求下,可得其迭代次数k=15,这说明改进后的Newton
法法速度确实比较快。
结论:
对于二分法,只要能够保证在给定的区间内有根,使能够收敛的,当时收敛的速度和给定
的区间有关,二且总体上来说速度比较慢。
Newton法,收敛速度要比二分法快,但是最终其
收敛的结果与初值的选取有关,初值不同,收敛的结果也可能不一样,也就是结果可能不时预
期需要得结果。
改进的Newton法求解重根问题时,如果初值不当,可能会不收敛,这一点非常重要,当然初值合适,相同情况下其速度要比Newton法快得多。
实验报告二
题目:
Gauss列主元消去法
摘要:
求解线性方程组的方法很多,主要分为直接法和间接法。
本实验运用直接法的Guass消
去法,并采用选主元的方法对方程组进行求解。
前言:
(目的和意义)
1.学习Gauss消去法的原理。
2.了解列主元的意义。
3.确定什么时候系数阵要选主元
数学原理:
由于一般线性方程在使用Gauss消去法求解时,从求解的过程中可以看到,若akk1)=O,
则必须进行行交换,才能使消去过程进行下去。
有的时候即使akk1)0,但是其绝对值非常小,由于机器舍入误差的影响,消去过程也会出现不稳定得现象,导致结果不正确。
因此有必要进行列主元技术,以最大可能的消除这种现象。
这一技术要寻找行r,使得
并将第r行和第k行的元素进行交换,以使得当前的akk1)的数值比0要大的多。
这种列主元的消去法的主要步骤如下:
1.消元过程
对k=1,2,…,n-1,进行如下步骤。
1)选主元,记
若|ark|很小,这说明方程的系数矩阵严重病态,给出警告,提示结果可能不对。
2)交换增广阵A的r,k两行的元素。
arjakj(j=k,…,n+1)
3)计算消元
aja0aikaq/akk(i=k+1,…,n;j=k+1,,n+1)
2.回代过程
对k=n,n-1,…,1,进行如下计算
至此,完成了整个方程组的求解。
程序设计:
本实验采用Matlab的M文件编写。
Gauss消去法源程序:
clear
a=input('输入系数阵:
>>\n')
b=input('输入列阵b:
>>\n')
n=length(b);
A=[ab]
x=zeros(n,1);
%%%函数主体
fork=1:
n-1;
%%%是否进行主兀选取
ifabs(A(k,k))vyipusilong;%事先给定的认为有必要选主元的小数yzhuyuan=1;
elseyzhuyuan=0;
end
ifyzhuyuan;
%%%%选主元
t=A(k,k);forr=k+1:
n;
ifabs(A(r,k))>abs(t)
p=r;
elsep=k;
end
end
%%%交换元素
ifp〜=k;
forq=k:
n+1;
s=A(k,q);
A(k,q)=A(p,q);A(p,q)=s;
end
end
end
%%%判断系数矩阵是否奇异或病态非常严重
ifabs(A(k,k))disp(矩阵奇异,解可能不正确’)
end
%%%%计算消元,得三角阵
forr=k+1:
n;
m=A(r,k)/A(k,k);
forq=k:
n+1;
A(r,q)=A(r,q)-A(k,q)*m;
end
end
end
%%%%求解x
x(n)=A(n,n+1)/A(n,n);
fork=n-1:
-1:
1;
s=0;
forr=k+1:
n;s=s+A(k,r)*x(r);
end
t=(A(k,n+1)-s)
x(k)=(A(k,n+1)-s)/A(k,k)
end
结果分析和讨论:
26x
例:
求解方程575y
321z
22
34。
其中为一小数,当
10
5101420
10,10,10,10时,分别
采用列主元和不列主元的Gauss消去法求解,并比较结果。
记Emax为求出的解代入方程后的最大误差,按要求,计算结果如下:
当105时,不选主元和选主元的计算结果如下,其中前一列为不选主元结果,后一列为选主元结果,下同。
Emax=,0
此时,由于不是很小,机器误差就不是很大,由Emax可以看出不选主元的计算结果精度还可以,因此此时可以考虑不选主元以减少计算量。
当1010时,不选主元和选主元的计算结果如下
Emax=,0
此时由Emax可以看出不选主元的计算精度就不好了,误差开始增大。
当1014时,不选主元和选主元的计算结果如下
0000000
Emax=,0
此时由Emax可以看出,不选主元的结果应该可以说是不正确了,这是由机器误差引起的。
当1020时,不选主元和选主元的计算结果如下
NaN1
NaN2
NaN3
Emax=NaN,0
不选主元时,程序报错:
Warning:
Dividebyzero.。
这是因为机器计算的最小精度为10-15,所以此时的1020就认为是0,故出现了错误现象。
而选主元时则没有这种现象,而且由Emax
可以看出选主元时的结果应该是精确解。
结论:
采用Gauss消去法时,如果在消元时对角线上的元素始终较大(假如大于10-5),那么本
方法不需要进行列主元计算,计算结果一般就可以达到要求,否则必须进行列主元这一步,以减少机器误差带来的影响,使方法得出的结果正确。
实验报告三
题目:
Rung现象产生和克服
摘要:
由于高次多项式插值不收敛,会产生Runge现象,本实验在给出具体的实例后,采用分段线性插值和三次样条插值的方法有效的克服了这一现象,而且还取的很好的插值效果。
前言:
(目的和意义)
1.深刻认识多项式插值的缺点。
2.明确插值的不收敛性怎样克服。
3.明确精度与节点和插值方法的关系。
数学原理:
在给定n+1个节点和相应的函数值以后构造n次的Lagrange插值多项式,实验结果表明
(见后面的图)这种多项式并不是随着次数的升高对函数的逼近越来越好,这种现象就是Rung
现象。
解决Rung现象的方法通常有分段线性插值、三次样条插值等方法。
分段线性插值:
设在区间[a,b]上,给定n+1个插值节点
a=xo和相应的函数值yo,yi,-•,yn,,求作一个插值函数(x),具有如下性质:
1)区)Yj,j=0,1,…,n。
2)(x)在每个区间[xi,刈上是线性连续函数。
则插值函数(x)称为区间[a,b]上对应n个
数据点的分段线性插值函数。
三次样条插值:
给定区间[a,b]一个分划
/:
a=X0若函数S(x)满足下列条件:
1)S(x)在每个区间[Xi,xj]上是不高于3次的多项式。
2)S(x)及其2阶导数在[a,b]上连续。
则称S(x)使关于分划/的三次样条函数。
程序设计:
本实验采用Matlab的M文件编写。
其中待插值的方程写成function的方式,如下
functiony=f(x);
y=1/(1+25*x*x);
写成如上形式即可,下面给出主程序
Lagrange插值源程序:
n=input('将区间分为的等份数输入:
\n');
s=[-1+2/n*[0:
n]];%%%给定的定点,Rf为给定的函数x=-1:
:
1;
f=0;
forq=1:
n+1;
l=1;%求插值基函数
fork=1:
n+1;
ifk〜=q;
l=l.*(x-s(k))./(s(q)-s(k));
else
1=1;
end
endf=f+Rf(s(q))*l;%求插值函数end
plot(x,f,'r')%作出插值函数曲线gridon
holdon
分段线性插值源程序
clear
n=input('将区间分为的等份数输入:
\n');s=[-1+2/n*[0:
n]];%%%给定的定点,Rf为给定的函数m=0;
hh=;
forx=-1:
hh:
1;
ff=0;
fork=1:
n+1;%%%求插值基函数
switchk
case1
ifx<=s
(2);l=(x-s
(2))./(s
(1)-s
(2));
else
l=0;
end
casen+1
ifx>s(n);
l=(x-s(n))./(s(n+1)-s(n));
else
l=0;
end
otherwise
ifx>=s(k-1)&x<=s(k);l=(x-s(k-1))./(s(k)-s(k-1));elseifx>=s(k)&x<=s(k+1);
l=(x-s(k+1))./(s(k)-s(k+1));
else
l=0;
end
end
end
ff=ff+Rf(s(k))*l;%%求插值函数值
end
m=m+1;f(m)=ff;
end
%%%作出曲线
x=-1:
hh:
1;
plot(x,f,'r');
gridon
holdon
三次样条插值源程序:
(采用第一边界条件)
clear
n=input('将区间分为的等份数输入:
\n');
%%%插值区间
a=-1;
b=1;
hh=;%画图的步长
s=[a+(b-a)/n*[0:
n]];%%%给定的定点,Rf为给定的函数%%%%第一边界条件Rf"(-1),Rf"
(1)v=5000*1/(1+25*a*a)A3-50/(1+25*a*a)A4;
fork=1:
n;%取出节点间距
h(k)=s(k+1)-s(k);
end
fork=1:
n-1;%求出系数向量lamuda,miula(k)=h(k+1)/(h(k+1)+h(k));
miu(k)=1-la(k);
end
%%%%赋值系数矩阵A
fork=1:
n-1;
forp=1:
n-1;
switchp
casek
A(k,p)=2;
casek1
A(k,p)=miu(p+1);
casek+1
A(k,p)=la(p-1);
otherwise
A(k,p)=0;
end
end
end
%%%%求出d阵
fork=1:
n-1;
switchk
case1
d(k)=6*f2c([s(k)s(k+1)s(k+2)])-miu(k)*v;
casen1
d(k)=6*f2c([s(k)s(k+1)s(k+2)])-la(k)*v;otherwise
d(k)=6*f2c([s(k)s(k+1)s(k+2)]);
end
end
%%%%求解M阵
M=A\d';
M=[v;M;v];
%%%%
m=0;
f=0;
forx=a:
hh:
b;
ifx==a;
p=1;
else
p=ceil((x-s
(1))/((b-a)/n));
end
ff仁0;
ff2=0;
ff3=0;
ff4=0;
m=m+1;
ff1=1/h(p)*(s(p+1)-x)A3*M(p)/6;
ff2=1/h(p)*(x-s(p))A3*M(p+1)/6;ff3=((Rf(s(p+1))-Rf(s(p)))/h(p)-h(p)*(M(p+1)-M(p))/6)*(x-s(p));ff4=Rf(s(p))-M(p)*h(p)*h(p)/6;
f(m)=ff1+ff2+ff3+ff4;
end
%%%作出插值图形
x=a:
hh:
b;
plot(x,f,'k')
holdon
gridon
结果分析和讨论:
本实验采用函数f(x)1一2进行数值插值,插值区间为卜1,1],给定节点为
125x
Xj=-1+jh,h=,j=0,…,n。
下面分别给出Lagrange插值,三次样条插值,线性插值的函数曲线和数据表。
图中只标出Lagrange插值的十次多项式的曲线,其它曲线没有标出,从数据表中可以看出具体的误差。
表中,Li0(x)为Lagrange插值的10次多项式,Sw(x),S40(x)分别代表n=10,40的三次样条插值函数,X10(x),X40(x)分别代表n=10,40的线性分段插值函数。
xf(x)L10(x)Sl0(x)S40(x)X10(x)X40(x)
0
从以上结果可以看到,用三次样条插值和线性分段插值,不会出现多项式插值是出现的Runge现象,插值效果明显提高。
进一步说,为了提高插值精度,用三次样条插值和线性分段插值是可以增加插值节点的办法来满足要求,而用多项式插值函数时,节点数的增加必然会使
多项式的次数增加,这样会引起数值不稳定,所以说这两种插值要比多项式插值好的多。
而且
在给定节点数的条件下,三次样条插值的精度要优于线性分段插值,曲线的光滑性也要好一些。
实验报告四
题目:
多项式最小二乘法
摘要:
对于具体实验时,通常不是先给出函数的解析式,再进行实验,而是通过实验的观察和测量给出离散的一些点,再来求出具体的函数解析式。
又因为测量误差的存在,实际真实的解
析式曲线并不一定通过测量给出的所有点。
最小二乘法是求解这一问题的很好的方法,本实验
运用这一方法实现对给定数据的拟合。
前言:
(目的和意义)
1.学习使用最小二成法的原理
2.了解法方程的特性数学原理:
对于给定的测量数据(xi,fi)(i=1,2,…,n),设函数分布为
特别的,取j(x)为多项式
j(x)xj(j=0,1,…,m)
则根据最小二乘法原理,可以构造泛函令
H
0(k=0,1,…,m)
ak
则可以得到法方程
求该解方程组,则可以得到解ao,a!
,am,因此可得到数据的最小二乘解
程序设计:
本实验采用Matlab的M文件编写。
其中多项式函数jxj写成function的方式,如下
functiony=fai(x,j)
y=1;
fori=1:
j
y=x.*y;
end
写成如上形式即可,下面给出主程序。
多项式最小二乘法源程序
clear
%%%给定测量数据点(s,f)s=[3456789];
f=[];
%%%计算给定的数据点的数目
n=length(f);
%%%给定需要拟合的数据的最高次多项式的次数
m=10;
%%%程序主体
fork=0:
m;
g=zeros(1,m+1);
forj=0:
m;
t=0;
fori=1:
n;%计算内积(fai(si),fai(si))t=t+fai(s(i),j)*fai(s(i),k);
endg(j+1)=t;
end
A(k+1,:
)=g;%法方程的系数矩阵
t=0;
fori=1:
n;%计算内积(f(si),fai(si))
t=t+f(i)*fai(s(i),k);
end
b(k+1,1)=t;
end
a=A\b%求出多项式系数
x=[s
(1):
:
s(n)]';
y=0;
fori=0:
m;
y=y+a(i+1)*fai(x,i);
end
plot(x,y)%作出拟合成的多项式的曲线
gridonholdon
plot(s,f,'rx')%在上图中标记给定的点
结果分析和讨论:
例用最小二乘法处理下面的实验数据.
Xi
3
4
5
6
7
8
9
fi
并作出f(x)的近似分布图。
分别采用一次,二次和五次多项式来拟合数据得到相应的拟合多项式为:
y1=+;
y2=+;
y5=+;
分别作出它们的曲线图,图中点划线为y1曲线,实线为y2曲线,虚线为y5曲线。
’x为给定的数据点。
从图中可以看出并不是多项式次数越高越好,次数高了,曲线越能给定点处和实际吻合,但别的地方就很差了。
因此,本例选用一次和两次的多项式拟合应该就可以了。
实验报告五
题目:
Romberg积分法
摘要:
对于实际的工程积分问题,很难应用NewtoeLeibnitz公式去求解。
因此应用数值方法进行求解积分问题已经有着很广泛的应用,本文基于Romberg积分法来解决一类积分问题。
前言:
(目的和意义)
1.理解和掌握Romberg积分法的原理;
2.学会使用Romberg积分法;
3.明确Romberg积分法的收敛速度及应用时容易出现的问题。
数学原理:
b
考虑积分1(f)f(x)dx,欲求其近似值,通常有复化的梯形公式、Simpsion公式和Cotes
a
公式。
但是给定一个精度,这些公式达到要求的速度很缓慢。
如何提高收敛速度,自然是人们极为关心的课题。
为此,记T1,k为将区间[a,b]进行2k等分的复化的梯形公式计算结果,记T2,k
为将区间[a,b]进行2k等分的复化的Simpsion公式计算结果,记T3,k为将区间[a,b]进行2k等分的复化的Cotes公式计算结果。
根据Richardson外推加速方法,可以得到收敛速度较快的
Romberg积分法。
其具体的计算公式为:
1.准备初值,计算
2.按梯形公式的递推关系,计算
3.按Romberg积分公式计算加速值
m1
Tm,k
m=2,…,k
4Tm1,k1mTm1,km
4.精度控制。
对给定的精度R,若
则终止计算,并取Tm,1为所求结果;否则返回2重复计算,直至满足要求的精度为止程序设计:
本实验采用Matlab的M文件编写。
其中待积分的函数写成function的方式,例如如下
functionyy=f(x,y);
yy=x.A3;
写成如上形式即可,下面给出主程序
Romberg积分法源程序%%%Romberg积分法clear
%%%积分区间
b=3;
a=1;
%%%精度要求
R=1e-5;