哈工大计算方法上机实验报告概要Word文档格式.docx
《哈工大计算方法上机实验报告概要Word文档格式.docx》由会员分享,可在线阅读,更多相关《哈工大计算方法上机实验报告概要Word文档格式.docx(32页珍藏版)》请在冰豆网上搜索。
0;
a=c;
else
b=c;
end
R=b-a;
%求出误差
k=k+1;
end
x=c%给出解
Newton法及改进的Newton法源程序:
%%%%输入函数
f=input('
请输入需要求解函数>
>
'
'
s'
)
%%%求解f(x)的导数
df=diff(f);
%%%改进常数或重根数
miu=2;
%%%初始值x0
x0=input('
inputinitialvaluex0>
);
%迭代次数
max=100;
%最大迭代次数
R=eval(subs(f,'
x0'
x'
));
%求解f(x0),以确定初值x0时否就是解
while(abs(R)>
1e-8)
x1=x0-miu*eval(subs(f,'
))/eval(subs(df,'
R=x1-x0;
x0=x1;
k=k+1;
if(eval(subs(f,'
))<
1e-10);
break
ifk>
max;
%如果迭代次数大于给定值,认为迭代不收敛,重新输入初值
ss=input('
mayberesultiserror,chooseanewx0,y/n?
ifstrcmp(ss,'
y'
x0=input('
k=0;
k;
%给出迭代次数
x=x0;
%给出解
结果分析和讨论:
1.用二分法计算方程
在[1,2]内的根。
(
下同)
计算结果为
x=1.40441513061523;
f(x)=-3.797205105904311e-007;
k=18;
由f(x)知结果满足要求,但迭代次数比较多,方法收敛速度比较慢。
2.用二分法计算方程
在[1,1.5]内的根。
x=1.32471847534180;
f(x)=2.209494846194815e-006;
k=17;
由f(x)知结果满足要求,但迭代次数还是比较多。
3.用Newton法求解下列方程
a)
x0=0.5;
x=0.56714329040978;
f(x)=2.220446049250313e-016;
k=4;
由f(x)知结果满足要求,而且又迭代次数只有4次看出收敛速度很快。
b)
x0=1;
c)
x0=0.45,x0=0.65;
当x0=0.45时,计算结果为
x=0.49999999999983;
f(x)=-8.362754932994584e-014;
由f(x)知结果满足要求,而且又迭代次数只有4次看出收敛速度很快,实际上该方程确实有真解x=0.5。
当x0=0.65时,计算结果为
x=0.50000000000000;
f(x)=0;
k=9;
由f(x)知结果满足要求,实际上该方程确实有真解x=0.5,但迭代次数增多,实际上当取x0〉0.68时,x≈1,就变成了方程的另一个解,这说明Newton法收敛与初值很有关系,有的时候甚至可能不收敛。
4.用改进的Newton法求解,有2重根,取
x0=0.55;
并与3.中的c)比较结果。
当x0=0.55时,程序死循环,无法计算,也就是说不收敛。
改
时,结果收敛为
x=0.50000087704286;
f(x)=4.385198907621127e-007;
k=16;
显然这个结果不是很好,而且也不是收敛至方程的2重根上。
当x0=0.85时,结果收敛为
x=1.00000000000489;
f(x)=2.394337647718737e-023;
这次达到了预期的结果,这说明初值的选取很重要,直接关系到方法的收敛性,实际上直接用Newton法,在给定同样的条件和精度要求下,可得其迭代次数k=15,这说明改进后的Newton法法速度确实比较快。
结论:
对于二分法,只要能够保证在给定的区间内有根,使能够收敛的,当时收敛的速度和给定的区间有关,二且总体上来说速度比较慢。
Newton法,收敛速度要比二分法快,但是最终其收敛的结果与初值的选取有关,初值不同,收敛的结果也可能不一样,也就是结果可能不时预期需要得结果。
改进的Newton法求解重根问题时,如果初值不当,可能会不收敛,这一点非常重要,当然初值合适,相同情况下其速度要比Newton法快得多。
实验报告二
Gauss列主元消去法
求解线性方程组的方法很多,主要分为直接法和间接法。
本实验运用直接法的Guass消去法,并采用选主元的方法对方程组进行求解。
1.学习Gauss消去法的原理。
2.了解列主元的意义。
3.确定什么时候系数阵要选主元
由于一般线性方程在使用Gauss消去法求解时,从求解的过程中可以看到,若
=0,则必须进行行交换,才能使消去过程进行下去。
有的时候即使
0,但是其绝对值非常小,由于机器舍入误差的影响,消去过程也会出现不稳定得现象,导致结果不正确。
因此有必要进行列主元技术,以最大可能的消除这种现象。
这一技术要寻找行r,使得
并将第r行和第k行的元素进行交换,以使得当前的
的数值比0要大的多。
这种列主元的消去法的主要步骤如下:
1.消元过程
对k=1,2,…,n-1,进行如下步骤。
1)选主元,记
若
很小,这说明方程的系数矩阵严重病态,给出警告,提示结果可能不对。
2)交换增广阵A的r,k两行的元素。
(j=k,…,n+1)
3)计算消元
(i=k+1,…,n;
j=k+1,……,n+1)
2.回代过程
对k=n,n-1,…,1,进行如下计算
至此,完成了整个方程组的求解。
Gauss消去法源程序:
a=input('
输入系数阵:
\n'
b=input('
输入列阵b:
n=length(b);
A=[ab]
x=zeros(n,1);
%%%函数主体
fork=1:
n-1;
%%%是否进行主元选取
ifabs(A(k,k))<
yipusilong;
%事先给定的认为有必要选主元的小数
yzhuyuan=1;
elseyzhuyuan=0;
ifyzhuyuan;
%%%%选主元
t=A(k,k);
forr=k+1:
n;
ifabs(A(r,k))>
abs(t)
p=r;
elsep=k;
%%%交换元素
ifp~=k;
forq=k:
n+1;
s=A(k,q);
A(k,q)=A(p,q);
A(p,q)=s;
end
%%%判断系数矩阵是否奇异或病态非常严重
yipusilong
disp(‘矩阵奇异,解可能不正确’)
%%%%计算消元,得三角阵
m=A(r,k)/A(k,k);
A(r,q)=A(r,q)-A(k,q)*m;
%%%%求解x
x(n)=A(n,n+1)/A(n,n);
fork=n-1:
-1:
1;
s=0;
s=s+A(k,r)*x(r);
t=(A(k,n+1)-s)
x(k)=(A(k,n+1)-s)/A(k,k)
例:
求解方程
。
其中
为一小数,当
时,分别采用列主元和不列主元的Gauss消去法求解,并比较结果。
记Emax为求出的解代入方程后的最大误差,按要求,计算结果如下:
当
时,不选主元和选主元的计算结果如下,其中前一列为不选主元结果,后一列为选主元结果,下同。
0.999999347683910.99999934782651
2.000002174219722.00000217391163
2.999997608594512.99999760869721
Emax=9.301857062382624e-010,0
此时,由于
不是很小,机器误差就不是很大,由Emax可以看出不选主元的计算结果精度还可以,因此此时可以考虑不选主元以减少计算量。
时,不选主元和选主元的计算结果如下
1.000017846308770.99999999999348
1.999980097208072.00000000002174
3.000006634247312.99999999997609
Emax=2.036758973744668e-005,0
此时由Emax可以看出不选主元的计算精度就不好了,误差开始增大。
1.421085471520201.00000000000000
1.666666666666662.00000000000000
3.11111111111111300000000000000
Emax=0.70770085900503,0
此时由Emax可以看出,不选主元的结果应该可以说是不正确了,这是由机器误差引起的。
NaN1
NaN2
NaN3
Emax=NaN,0
不选主元时,程序报错:
Warning:
Dividebyzero.。
这是因为机器计算的最小精度为10-15,所以此时的
就认为是0,故出现了错误现象。
而选主元时则没有这种现象,而且由Emax可以看出选主元时的结果应该是精确解。
采用Gauss消去法时,如果在消元时对角线上的元素始终较大(假如大于10-5),那么本方法不需要进行列主元计算,计算结果一般就可以达到要求,否则必须进行列主元这一步,以减少机器误差带来的影响,使方法得出的结果正确。
实验报告三
Rung现象产生和克服
由于高次多项式插值不收敛,会产生Runge现象,本实验在给出具体的实例后,采用分段线性插值和三次样条插值的方法有效的克服了这一现象,而且还取的很好的插值效果。
1.深刻认识多项式插值的缺点。
2.明确插值的不收敛性怎样克服。
3.明确精度与节点和插值方法的关系。
在给定n+1个节点和相应的函数值以后构造n次的Lagrange插值多项式,实验结果表明(见后面的图)这种多项式并不是随着次数的升高对函数的逼近越来越好,这种现象就是Rung现象。
解决Rung现象的方法通常有分段线性插值、三次样条插值等方法。
分段线性插值:
设在区间[a,b]上,给定n+1个插值节点
a=x0<
x1<
…<
xn=b
和相应的函数值y0,y1,…,yn,,求作一个插值函数
,具有如下性质:
1)
,j=0,1,…,n。
2)
在每个区间[xi,xj]上是线性连续函数。
则插值函数
称为区间[a,b]上对应n个数据点的分段线性插值函数。
三次样条插值:
给定区间[a,b]一个分划
⊿:
xN=b
若函数S(x)满足下列条件:
1)S(x)在每个区间[xi,xj]上是不高于3次的多项式。
2)S(x)及其2阶导数在[a,b]上连续。
则称S(x)使关于分划⊿的三次样条函数。
其中待插值的方程写成function的方式,如下
y=1/(1+25*x*x);
写成如上形式即可,下面给出主程序
Lagrange插值源程序:
n=input('
将区间分为的等份数输入:
s=[-1+2/n*[0:
n]];
%%%给定的定点,Rf为给定的函数
x=-1:
0.01:
f=0;
forq=1:
l=1;
%求插值基函数
fork=1:
ifk~=q;
l=l.*(x-s(k))./(s(q)-s(k));
l=l;
f=f+Rf(s(q))*l;
%求插值函数
plot(x,f,'
r'
)%作出插值函数曲线
gridon
holdon
分段线性插值源程序
m=0;
hh=0.001;
forx=-1:
hh:
ff=0;
%%%求插值基函数
switchk
case1
ifx<
=s
(2);
l=(x-s
(2))./(s
(1)-s
(2));
l=0;
casen+1
ifx>
s(n);
l=(x-s(n))./(s(n+1)-s(n));
otherwise
=s(k-1)&
x<
=s(k);
l=(x-s(k-1))./(s(k)-s(k-1));
elseifx>
=s(k)&
=s(k+1);
l=(x-s(k+1))./(s(k)-s(k+1));
else
ff=ff+Rf(s(k))*l;
%%求插值函数值
m=m+1;
f(m)=ff;
end
%%%作出曲线
gridon
holdon
三次样条插值源程序:
(采用第一边界条件)
%%%插值区间
a=-1;
b=1;
%画图的步长
s=[a+(b-a)/n*[0:
%%%%第一边界条件Rf"
(-1),Rf"
(1)
v=5000*1/(1+25*a*a)^3-50/(1+25*a*a)^4;
%取出节点间距
h(k)=s(k+1)-s(k);
%求出系数向量lamuda,miu
la(k)=h(k+1)/(h(k+1)+h(k));
miu(k)=1-la(k);
%%%%赋值系数矩阵A
forp=1:
switchp
casek
A(k,p)=2;
casek-1
A(k,p)=miu(p+1);
casek+1
A(k,p)=la(p-1);
otherwise
A(k,p)=0;
%%%%求出d阵
d(k)=6*f2c([s(k)s(k+1)s(k+2)])-miu(k)*v;
casen-1
d(k)=6*f2c([s(k)s(k+1)s(k+2)])-la(k)*v;
d(k)=6*f2c([s(k)s(k+1)s(k+2)]);
%%%%求解M阵
M=A\d'
;
M=[v;
M;
v];
%%%%
forx=a:
b;
ifx==a;
p=1;
p=ceil((x-s
(1))/((b-a)/n));
ff1=0;
ff2=0;
ff3=0;
ff4=0;
ff1=1/h(p)*(s(p+1)-x)^3*M(p)/6;
ff2=1/h(p)*(x-s(p))^3*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;
%%%作出插值图形
x=a:
k'
本实验采用函数
进行数值插值,插值区间为[-1,1],给定节点为
xj=-1+jh,h=0.1,j=0,…,n。
下面分别给出Lagrange插值,三次样条插值,线性插值的函数曲线和数据表。
图中只标出Lagrange插值的十次多项式的曲线,其它曲线没有标出,从数据表中可以看出具体的误差。
表中,L10(x)为Lagrange插值的10次多项式,S10(x),S40(x)分别代表n=10,40的三次样条插值函数,X10(x),X40(x)分别代表n=10,40的线性分段插值函数。
xf(x)L10(x)S10(x)S40(x)X10(x)X40(x)
-1.000000000000000.038461538461540.038461538461540.038461538461540.038461538461540.038461538461540.03846153846154
-0.950000000000000.042440318302391.923631149719200.042408331510400.042440318302390.043552036199100.04244031830239
-0.900000000000000.047058823529411.578720990349260.047096975854580.047058823529410.048642533936650.04705882352941
-0.850000000000000.052459016393440.719459128379820.052558399239790.052459016393440.053733031674210.05245901639344
-0.800000000000000.058823529411760.058823529411760.058823529411760.058823529411760.058823529411760.05882352941176
-0.750000000000000.06639004149378-0.231461749896740.066039861727440.066390041493780.069117647058820.06639004149378
-0.700000000000000.07547169811321-0.226196289062500.074821161988660.075471698113210.079411764705880.07547169811321
-0.650000000000000.08648648648649-0.072604203224180.085897763608490.086486486486490.089705882352940.08648648648649
-0.600000000000000.100000000000000.100000000000000.100000000000000.100000000000000.100000000000000.10000000000000
-0.550000000000000.116788321167880.215591878912570.117838330177130.116788321167880.125000000000000.11678832116788
-0.500000000000000.137********2760.253755457261030.140043715557300.137********2760.150********0000.137********276
-0.450000000000000.164948453608250.234968543052670.167227243158830.164948453608250.175********0000.16494845360825
-0.400000000000000.200000000000000.200000000000000.200000000000000.200000000000000.200000000000000.20000000000000
-0.350000000000000.246153846153850.190580466753760.240547994034640.246153846153850.275000000000000.24615384615385
-0.300000000000000.307692307692310.235346591310800.297356916958600.307692307692310.3500