哈工大计算方法上机实验报告1.docx
- 文档编号:2815511
- 上传时间:2022-11-15
- 格式:DOCX
- 页数:32
- 大小:186.96KB
哈工大计算方法上机实验报告1.docx
《哈工大计算方法上机实验报告1.docx》由会员分享,可在线阅读,更多相关《哈工大计算方法上机实验报告1.docx(32页珍藏版)》请在冰豆网上搜索。
哈工大计算方法上机实验报告1
实验报告一
题目:
非线性方程求解
摘要:
非线性方程的解析解通常很难给出,因此线性方程的数值解法就尤为重要。
本实验采用两种常见的求解方法二分法和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法通常预先要给出一个猜测初值x0,然后根据其迭代公式
产生逼近解x*的迭代数列{xk},这就是Newton法的思想。
当x0接近x*时收敛很快,但是当x0选择不好时,可能会发散,因此初值的选取很重要。
另外,若将该迭代公式改进为
其中r为要求的方程的根的重数,这就是改进的Newton法,当求解已知重数的方程的根时,在同种条件下其收敛速度要比Newton法快的多。
程序设计:
本实验采用Matlab的M文件编写。
其中待求解的方程写成function的方式,如下
functiony=f(x);
y=-x*x-sin(x);
写成如上形式即可,下面给出主程序。
二分法源程序:
clear
%%%给定求解区间
b=1.5;
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
x0=input('inputinitialvaluex0>>');
k=0;%迭代次数
max=100;%最大迭代次数
R=eval(subs(f,'x0','x'));%求解f(x0),以确定初值x0时否就是解
while(abs(R)>1e-8)
x1=x0-miu*eval(subs(f,'x0','x'))/eval(subs(df,'x0','x'));
R=x1-x0;
x0=x1;
k=k+1;
if(eval(subs(f,'x0','x'))<1e-10);
break
end
ifk>max;%如果迭代次数大于给定值,认为迭代不收敛,重新输入初值
ss=input('mayberesultiserror,chooseanewx0,y/n?
>>','s');
ifstrcmp(ss,'y')
x0=input('inputinitialvaluex0>>');
k=0;
else
break
end
end
end
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;
k=4;
由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;
k=4;
这次达到了预期的结果,这说明初值的选取很重要,直接关系到方法的收敛性,实际上直接用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,进行如下计算
至此,完成了整个方程组的求解。
程序设计:
本实验采用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)) 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 结果分析和讨论: 例: 求解方程 。 其中 为一小数,当 时,分别采用列主元和不列主元的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,所以
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 哈工大 计算方法 上机 实验 报告