(完整版)哈工大-数值分析上机实验报告
数值分析上机实验报告3
实验报告三题目:函数逼近——曲线拟合目的:掌握曲线拟合基本使用方法数学原理:[P,S]=polyfit(x,y,3)其中x,y为取样值,3为得出的结果的最高次数。
P为对应次数的系数,S为误差值向量,其中x,y是等长的向量,P是一个长度为m+1的向量。
结果分析和讨论:23.观察物体的运动,得出时间t与距离s的关系如表,求运动方程。
t=[0,0.9,1.9,3.0,3.9,5.0];s=[0,10,30,50,80,110];[P,S]=polyfit(t,s,5)P =-0.5432 6.4647 -26.5609 46.1436 -13.2601 -0.0000S =R: [6x6 double]df: 0normr: 1.2579e-012所以得到方程为:-13.2601x46.1436x-26.5609x6.4647x-0.5432x2345++24.在某化学反应堆里,根据实验所得分解物的质量分数y与时间t的关系,用最小拟合求y=F(t);>> x=0:5:55;y=[0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.62,4.64];>> [P,S]=polyfit(x,y,5)P =0.0000 -0.0000 0.0002 -0.0084 0.2851 0.0082S =R: [6x6 double]df: 6normr: 0.0487所以得到方程为:0082.02851.00084.00002.023++-xxx结论:在23题中计算的结果误差为4.5769,而在24中计算的结果误差为0.0487,说明对于曲线拟合来说,总会有误差,因为取样点并不是都过拟合的曲线的。
数值分析上机实习报告
指导教师:姓名:学号:专业:联系电话:上海交通大学目录序言 (3)实验课题(一) 雅可比迭代法和高斯-塞得尔迭代法的收敛性和收敛速度 (4)数值分析 (6)实验课题(二) 松弛因子对SOR法收敛速度的影响 (6)数值分析 (12)总结 (13)附录(程序清单) (14)1.雅可比迭代法和高斯-塞得尔迭代法的收敛性和收敛速度 (14)雅可比迭代法: (14)高斯-塞得尔迭代法: (16)2.松弛因子对SOR法收敛速度的影响 (18)松弛法(SOR) (18)序言随着科学技术的发展,提出了大量复杂的数值计算问题,在实际解决这些计算问题的长期过程中,形成了计算方法这门学科,专门研究各种数学问题的数值解法(近似解法),包括方法的构造和求解过程的误差分析,是一门内容丰富,有自身理论体系的实用性很强的学科。
解决工程问题,往往需要处理很多数学模型,这就要花费大量的人力和时间,但是还有不少数学模型无法用解析法得到解。
使用数值方法并利用计算机,就可以克服这些困难。
事实上,科学计算已经与理论分析、科学实验成为平行的研究和解决科技问题的科学手段,经常被科技工作者所采用。
作为科学计算的核心内容——数值分析(数值计算方法),已逐渐成为广大科技工作者必备的基本知识并越来越被人重视。
由于数值方法是解数值问题的系列计算公式,所以数值方法是否有效,不但与方法本身的好坏有关,而且与数值问题本身的好坏也有关,因此,研究数值方法时,不但需要研究数值方法的好坏,即数值稳定性问题,而且还需要研究数值问题本身的好坏,即数值问题的性态,以及它们的判别问题。
数值计算的绝大部分方法都具有近似性,而其理论又具有严密的科学性,方法的近似值正是建立在理论的严密性基础上,根据计算方法的这一特点。
因此不仅要求掌握和使用算法,还要重视必要的误差分析,以保证计算结果的可靠性。
数值计算还具有应用性强的特点,计算方法的绝大部分方法如求微分方程近似解,求积分近似值,求解超越方程,解线性方程组等都具有较强的实用性,而插值法,最小二乘法,样条函数等也都是工程技术领域中常用的,有实际应用价值的方法。
优质文档精选——数值分析上机实验报告
数值分析上机实验报告《数值分析》上机实验报告1.用Newton 法求方程 X 7-X 4+14=0在(0.1,1.9)中的近似根(初始近似值取为区间端点,迭代6次或误差小于0.00001)。
1.1 理论依据:设函数在有限区间[a ,b]上二阶导数存在,且满足条件{}αϕ上的惟一解在区间平方收敛于方程所生的迭代序列迭代过程由则对任意初始近似值达到的一个中使是其中上不变号在区间],[0)(3,2,1,0,)(')()(],,[x |))(),((|,|,)(||)(|.4;0)(.3],[)(.20)()(.110......b a x f x k x f x f x x x Newton b a b f a f mir b a c x f ab c f x f b a x f b f x f k k k k k k ==-==∈≤-≠>+令)9.1()9.1(0)8(4233642)(0)16(71127)(0)9.1(,0)1.0(,1428)(3225333647>⋅''<-=-=''<-=-='<>+-=f f x x x x x f x x x x x f f f x x x f故以1.9为起点⎪⎩⎪⎨⎧='-=+9.1)()(01x x f x f x x k k k k 如此一次一次的迭代,逼近x 的真实根。
当前后两个的差<=ε时,就认为求出了近似的根。
本程序用Newton 法求代数方程(最高次数不大于10)在(a,b )区间的根。
1.2 C 语言程序原代码:#include<stdio.h>#include<math.h> main(){double x2,f,f1;double x1=1.9; //取初值为 1.9 do{x2=x1;f=pow(x2,7)-28*pow(x2,4)+14; f1=7*pow(x2,6)-4*28*pow(x2,3); x1=x2-f/f1;}while(fabs(x1-x2)>=0.00001||x1<0.1); //限制循环次数 printf("计算结果:x=%f\n",x1);}1.3 运行结果:1.4 MATLAB 上机程序function y=Newton(f,df,x0,eps,M) d=0;for k=1:Mif feval(df,x0)==0d=2;breakelsex1=x0-feval(f,x0)/feval(df,x0);ende=abs(x1-x0);x0=x1;if e<=eps&&abs(feval(f,x1))<=epsd=1;breakendendif d==1y=x1;elseif d==0y='迭代M次失败';elsey= '奇异'endfunction y=df(x)y=7*x^6-28*4*x^3;Endfunction y=f(x)y=x^7-28*x^4+14;End>> x0=1.9;>> eps=0.00001;>> M=100;>> x=Newton('f','df',x0,eps,M);>> vpa(x,7)1.5 问题讨论:1.使用此方法求方解,用误差来控制循环迭代次数,可以在误差允许的范围内得到比较理想的计算结果。
数值分析上机实验报告(插值)
数值分析第一次上机练习实验报告——Lagrange 插值与三次样条插值一、 问题的描述设()2119f x x =+, []1,1x ∈-,取15iix =-+,0,1,2,...,10i =.试求出10次Lagrange 插值多项式()10L x 和三次样条插值函数()S x (采用自然边界条件),并用图画出()f x ,()10L x ,()S x .二、 方法描述——Lagrange 插值与三次样条插值我们取15i ix =-+,0,1,2,...,10i =,通过在i x 点的函数值()2119i i f x x =+来对原函数进行插值,我们记插值函数为()g x ,要求它满足如下条件:()()21,0,1,2,...,1019i i i g x f x i x ===+ (1)我们在此处要分别通过Lagrange 插值(即多项式插值)与三次样条插值的方法对原函数()2119f x x=+进行插值,看两种方法的插值结果,并进行结果的比较。
10次的Lagrange 插值多项式为:()()10100i i i L x y l x ==∑ (2)其中:()21,0,1,2,...,1019i i iy f x i x ===+ 以及()()()()()()()()()011011......,0,1,2,...,10......i i n i i i i i i i n x x x x x x x x l x i x x x x x x x x -+-+----==----我们根据(2)进行程序的编写,我们可以通过几个循环很容易实现函数的Lagrange 插值。
理论上我们根据区间[]1,1-上给出的节点做出的插值多项式()n L x 近似于()f x ,而多项式()n L x 的次数n 越高逼近()f x 的精度就越好。
但实际上并非如此,而是对任意的插值节点,当n →+∞的时候()n L x 不一定收敛到()f x ;而是有时会在插值区间的两端点附近会出现严重的()n L x 偏离()f x 的现象,即所谓的Runge 现象。
哈工大数值分析报告上机实验报告
实验报告一题目: Gauss 列主元消去法摘要:求解线性方程组地方法很多,主要分为直接法和间接法.本实验运用直接法地Guass 消去法,并采用选主元地方法对方程组进行求解.前言:(目地和意义)1. 学习Gauss 消去法地原理.2. 了解列主元地意义.3. 确定什么时候系数阵要选主元数学原理:由于一般线性方程在使用Gauss 消去法求解时,从求解地过程中可以看到,若)1(-k kk a =0,则必须进行行交换,才能使消去过程进行下去.有地时候即使≠-)1(k kk a 0,但是其绝对值非常小,由于机器舍入误差地影响,消去过程也会出现不稳定得现象,导致结果不正确.因此有必要进行列主元技术,以最大可能地消除这种现象.这一技术要寻找行r ,使得)1()1(max ||->-=k ik ki k rk a a 并将第r 行和第k 行地元素进行交换,以使得当前地)1(-k kk a 地数值比0要大地多.这种列主元地消去法地主要步骤如下:1. 消元过程对k =1,2,…,n -1,进行如下步骤.1) 选主元,记ik ki rk a a >=max || 若||rk a 很小,这说明方程地系数矩阵严重病态,给出警告,提示结果可能不对.2) 交换增广阵A 地r ,k 两行地元素.kj rj a a ↔ (j=k,…,n +1)3) 计算消元kk kj ik ij ij a a a a a /-= (i=k+1,…,n ; j =k +1,……,n +1)2. 回代过程对k = n , n -1,…,1,进行如下计算)/(11,∑-=+-=nk j kk j kj n k k a x a a x至此,完成了整个方程组地求解.程序设计:本实验采用Matlab地M文件编写.Gauss消去法源程序:cleara=input('输入系数阵:>>\n')b=input('输入列阵b:>>\n')n=length(b);A=[a b]x=zeros(n,1);%%%函数主体for k=1:n-1;%%%是否进行主元选取if abs(A(k,k))<yipusilong;%事先给定地认为有必要选主元地小数yzhuyuan=1;else yzhuyuan=0;endif yzhuyuan;%%%%选主元t=A(k,k);for r=k+1:n;if abs(A(r,k))>abs(t)p=r;else p=k;endend%%%交换元素if p~=k;for q=k:n+1;s=A(k,q);A(k,q)=A(p,q);A(p,q)=s;endendend%%%判断系数矩阵是否奇异或病态非常严重if abs(A(k,k))< yipusilongdisp(‘矩阵奇异,解可能不正确’)end%%%%计算消元,得三角阵for r=k+1:n;m=A(r,k)/A(k,k);for q=k:n+1;A(r,q)=A(r,q)-A(k,q)*m;endendend%%%%求解xx(n)=A(n,n+1)/A(n,n);for k=n-1:-1:1;s=0;for r=k+1:n;s=s+A(k,r)*x(r);endt=(A(k,n+1)-s)x(k)=(A(k,n+1)-s)/A(k,k)end结果分析和讨论:例:求解方程⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡321643.5072.12-623.4712.31-32108-z y x .求解地结果为:x =[]367257386.0,05088607.0-49105822.0-, 例:求解方程⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡73109104-10172-42-4z y x 求得地结果为:x =[]857142857.1,89285714.0-196428571.0, 结论:采用Gauss 消去法时,如果在消元时对角线上地元素始终较大(假如大于10-5),那么本方法不需要进行列主元计算,计算结果一般就可以达到要求,否则必须进行列主元这一步,以减少机器误差带来地影响,使方法得出地结果正确.实验报告二题目: Rung 现象产生和克服摘要:由于高次多项式插值不收敛,会产生Runge 现象,本实验在给出具体地实例后,采用分段线性插值和三次样条插值地方法有效地克服了这一现象,而且还取地很好地插值效果.前言:(目地和意义)1. 深刻认识多项式插值地缺点.2. 明确插值地不收敛性怎样克服.3. 明确精度与节点和插值方法地关系.数学原理:在给定n+1个节点和相应地函数值以后构造n 次地Lagrange 插值多项式,实验结果表明(见后面地图)这种多项式并不是随着次数地升高对函数地逼近越来越好,这种现象就是Rung 现象.解决Rung 现象地方法通常有分段线性插值、三次样条插值等方法.分段线性插值:设在区间[a, b ]上,给定n+1个插值节点a=x 0<x 1<…<x n =b和相应地函数值y 0,y 1,…,y n ,,求作一个插值函数)(x φ,具有如下性质:1) j j y x =)(φ,j=0,1,…,n .2) )(x φ在每个区间[x i , x j ]上是线性连续函数.则插值函数)(x φ称为区间[a, b ]上对应n 个数据点地分段线性插值函数.三次样条插值:给定区间[a, b ]一个分划⊿:a=x 0<x 1<…<x N =b若函数S(x)满足下列条件:1) S(x)在每个区间[x i , x j ]上是不高于3次地多项式.2) S(x)及其2阶导数在[a, b ]上连续.则称S(x)使关于分划⊿地三次样条函数. 程序设计流程:本实验采用Matlab 地M 文件编写.其中待插值地方程写成function 地方式,如下function y=f(x);y=1/(1+25*x*x );写成如上形式即可,下面给出主程序Lagrange 插值源程序:n=input('将区间分为地等份数输入:\n');s=[-1+2/n*[0:n]];%%%给定地定点,Rf为给定地函数x=-1:0.01:1;f=0;for q=1:n+1;l=1;%求插值基函数for k=1:n+1;if k~=q;l=l.*(x-s(k))./(s(q)-s(k));elsel=l;endendf=f+Rf(s(q))*l;%求插值函数endplot(x,f,'r')%作出插值函数曲线grid onhold on分段线性插值源程序clearn=input('将区间分为地等份数输入:\n');s=[-1+2/n*[0:n]];%%%给定地定点,Rf为给定地函数m=0;hh=0.001;for x=-1:hh:1;ff=0;for k=1:n+1;%%%求插值基函数switch kcase 1if x<=s(2);l=(x-s(2))./(s(1)-s(2));elsel=0;endcase n+1if x>s(n);l=(x-s(n))./(s(n+1)-s(n));elsel=0;endotherwiseif x>=s(k-1)&x<=s(k);l=(x-s(k-1))./(s(k)-s(k-1));else if x>=s(k)&x<=s(k+1);l=(x-s(k+1))./(s(k)-s(k+1));elsel=0;endendendff=ff+Rf(s(k))*l;%%求插值函数值endm=m+1;f(m)=ff;end%%%作出曲线x=-1:hh:1;plot(x,f,'r');grid onhold on三次样条插值源程序:(采用第一边界条件)clearn=input('将区间分为地等份数输入:\n');%%%插值区间a=-1;b=1;hh=0.001;%画图地步长s=[a+(b-a)/n*[0:n]];%%%给定地定点,Rf为给定地函数%%%%第一边界条件Rf"(-1),Rf"(1)v=5000*1/(1+25*a*a)^3-50/(1+25*a*a)^4;for k=1:n;%取出节点间距h(k)=s(k+1)-s(k);endfor k=1:n-1;%求出系数向量lamuda,miula(k)=h(k+1)/(h(k+1)+h(k));miu(k)=1-la(k);end%%%%赋值系数矩阵Afor k=1:n-1;for p=1:n-1;switch pcase kA(k,p)=2;case k-1A(k,p)=miu(p+1);case k+1A(k,p)=la(p-1);otherwiseA(k,p)=0;endendend%%%%求出d阵for k=1:n-1;switch kcase 1d(k)=6*f2c([s(k) s(k+1) s(k+2)])-miu(k)*v;case n-1d(k)=6*f2c([s(k) s(k+1) s(k+2)])-la(k)*v;otherwised(k)=6*f2c([s(k) s(k+1) s(k+2)]);endend%%%%求解M阵M=A\d';M=[v;M;v];%%%%m=0;f=0;for x=a:hh:b;if x==a;p=1;elsep=ceil((x-s(1))/((b-a)/n));endff1=0;ff2=0;ff3=0;ff4=0;m=m+1;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 ;end%%%作出插值图形x=a:hh:b;plot(x,f,'k')hold ongrid on结果分析和讨论: 本实验采用函数22511)(xx f +=进行数值插值,插值区间为[-1,1],给定节点为 x j =-1+jh ,h=0.1,j =0,…,n .下面分别给出Lagrang e 插值,三次样条插值,线性插值地函数曲线和数据表.图中只标出Lagrang e 插值地十次多项式地曲线,其它曲线没有标出,从数据表中可以看出具体地误差.表中,L10(x)为Lagrang e插值地10次多项式,S10(x),S40(x)分别代表n=10,40地三次样条插值函数,X10(x),X40(x)分别代表n=10,40地线性分段插值函数.x f(x)L10(x)S10(x) S40(x) X10(x) X40(x) -1.00000000000000 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154 -0.95000000000000 0.04244031830239 1.92363114971920 0.04240833151040 0.04244031830239 0.04355203619910 0.04244031830239 -0.90000000000000 0.04705882352941 1.57872099034926 0.04709697585458 0.04705882352941 0.04864253393665 0.04705882352941 -0.85000000000000 0.05245901639344 0.71945912837982 0.05255839923979 0.05245901639344 0.05373303167421 0.05245901639344 -0.80000000000000 0.05882352941176 0.05882352941176 0.05882352941176 0.05882352941176 0.05882352941176 0.05882352941176 -0.75000000000000 0.06639004149378 -0.23146174989674 0.06603986172744 0.06639004149378 0.06911764705882 0.06639004149378 -0.70000000000000 0.07547169811321 -0.22619628906250 0.07482116198866 0.07547169811321 0.07941176470588 0.07547169811321 -0.65000000000000 0.08648648648649 -0.07260420322418 0.08589776360849 0.08648648648649 0.08970588235294 0.08648648648649 -0.60000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.10000000000000 -0.55000000000000 0.11678832116788 0.21559187891257 0.11783833017713 0.11678832116788 0.12500000000000 0.11678832116788 -0.50000000000000 0.13793103448276 0.25375545726103 0.14004371555730 0.13793103448276 0.15000000000000 0.13793103448276 -0.45000000000000 0.16494845360825 0.23496854305267 0.16722724315883 0.16494845360825 0.17500000000000 0.16494845360825 -0.40000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.20000000000000 -0.35000000000000 0.24615384615385 0.19058046675376 0.24054799403464 0.24615384615385 0.27500000000000 0.24615384615385 -0.30000000000000 0.30769230769231 0.23534659131080 0.29735691695860 0.30769230769231 0.35000000000000 0.30769230769231 -0.25000000000000 0.39024390243902 0.34264123439789 0.38048738140327 0.39024390243902 0.42500000000000 0.39024390243902 -0.20000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.50000000000000 -0.15000000000000 0.64000000000000 0.67898957729340 0.65746969368431 0.64000000000000 0.62500000000000 0.64000000000000 -0.10000000000000 0.80000000000000 0.84340742982890 0.82052861660828 0.80000000000000 0.75000000000000 0.80000000000000 -0.05000000000000 0.94117647058824 0.95862704866073 0.94832323122810 0.94117647058824 0.87500000000000 0.941176470588240 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.000000000000001.000000000000000.05000000000000 0.94117647058824 0.95862704866073 0.94832323122810 0.94117647058824 0.87500000000000 0.941176470588240.10000000000000 0.80000000000000 0.84340742982890 0.82052861660828 0.80000000000000 0.75000000000000 0.800000000000000.15000000000000 0.64000000000000 0.67898957729340 0.65746969368431 0.64000000000000 0.62500000000000 0.640000000000000.20000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.50000000000000 0.500000000000000.25000000000000 0.39024390243902 0.34264123439789 0.38048738140327 0.39024390243902 0.42500000000000 0.390243902439020.30000000000000 0.30769230769231 0.23534659131080 0.29735691695860 0.30769230769231 0.350000000000000.307692307692310.35000000000000 0.24615384615385 0.19058046675376 0.24054799403464 0.24615384615385 0.27500000000000 0.246153846153850.40000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.20000000000000 0.200000000000000.45000000000000 0.16494845360825 0.23496854305267 0.16722724315883 0.16494845360825 0.17500000000000 0.164948453608250.50000000000000 0.13793103448276 0.25375545726103 0.14004371555730 0.13793103448276 0.15000000000000 0.137931034482760.55000000000000 0.11678832116788 0.21559187891257 0.11783833017713 0.11678832116788 0.12500000000000 0.116788321167880.60000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.10000000000000 0.100000000000000.65000000000000 0.08648648648649 -0.07260420322418 0.08589776360849 0.08648648648649 0.08970588235294 0.086486486486490.70000000000000 0.07547169811321 -0.22619628906250 0.07482116198866 0.07547169811321 0.07941176470588 0.075471698113210.75000000000000 0.06639004149378 -0.23146174989674 0.06603986172744 0.06639004149378 0.06911764705882 0.066390041493780.80000000000000 0.05882352941176 0.05882352941176 0.05882352941176 0.05882352941176 0.05882352941176 0.058823529411760.85000000000000 0.05245901639344 0.71945912837982 0.05255839923979 0.05245901639344 0.05373303167421 0.052459016393440.90000000000000 0.04705882352941 1.57872099034926 0.04709697585458 0.04705882352941 0.04864253393665 0.047058823529410.95000000000000 0.04244031830239 1.92363114971920 0.04240833151040 0.04244031830239 0.04355203619910 0.042440318302391.00000000000000 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154 0.03846153846154从以上结果可以看到,用三次样条插值和线性分段插值,不会出现多项式插值是出现地Runge现象,插值效果明显提高.进一步说,为了提高插值精度,用三次样条插值和线性分段插值是可以增加插值节点地办法来满足要求,而用多项式插值函数时,节点数地增加必然会使多项式地次数增加,这样会引起数值不稳定,所以说这两种插值要比多项式插值好地多.而且在给定节点数地条件下,三次样条插值地精度要优于线性分段插值,曲线地光滑性也要好一些.实验报告三题目: 多项式最小二乘法摘要:对于具体实验时,通常不是先给出函数地解析式,再进行实验,而是通过实验地观察和测量给出离散地一些点,再来求出具体地函数解析式.又因为测量误差地存在,实际真实地解析式曲线并不一定通过测量给出地所有点.最小二乘法是求解这一问题地很好地方法,本实验运用这一方法实现对给定数据地拟合. 前言:(目地和意义)1. 学习使用最小二成法地原理2. 了解法方程地特性 数学原理:对于给定地测量数据(x i ,f i )(i=1,2,…,n ),设函数分布为∑==mj j j x a x y 0)()(ϕ特别地,取)(x j ϕ为多项式j j x x =)(ϕ (j=0, 1,…,m )则根据最小二乘法原理,可以构造泛函∑∑==-=n i mj i j j i m x a f a a a H 110))((),,,(ϕ令0=∂∂ka H(k=0, 1,…,m ) 则可以得到法方程⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡),(),(),(),(),(),(),(),(),(),(),(),(1010101111000100m m m m m m m m f f f a a a ϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕ求该解方程组,则可以得到解m a a a ,,,10 ,因此可得到数据地最小二乘解∑=≈mj j j x a x f 0)()(ϕ程序设计:本实验采用Matlab 地M 文件编写.其中多项式函数j j x =ϕ写成function 地方式,如下function y=fai(x,j)y=1;for i=1:jy=x.*y;end写成如上形式即可,下面给出主程序.多项式最小二乘法源程序clear%%%给定测量数据点(s,f)s=[3 4 5 6 7 8 9];f=[2.01 2.98 3.50 5.02 5.47 6.02 7.05];%%%计算给定地数据点地数目n=length(f);%%%给定需要拟合地数据地最高次多项式地次数m=10;%%%程序主体for k=0:m;g=zeros(1,m+1);for j=0:m;t=0;for i=1:n;%计算内积(fai(si),fai(si))t=t+fai(s(i),j)*fai(s(i),k);endg(j+1)=t;endA(k+1,:)=g;%法方程地系数矩阵t=0;for i=1:n;%计算内积(f(si),fai(si))t=t+f(i)*fai(s(i),k);endb(k+1,1)=t;enda=A\b%求出多项式系数x=[s(1):0.01:s(n)]';y=0;for i=0:m;y=y+a(i+1)*fai(x,i);endplot(x,y)%作出拟合成地多项式地曲线grid onhold onplot(s,f,'rx') %在上图中标记给定地点结果分析和讨论:例用最小二乘法处理下面地实验数据.并作出)f地近似分布图.(x分别采用一次,二次和五次多项式来拟合数据得到相应地拟合多项式为:y1=-0.38643+0.82750x;y2=-1.03024+1.06893x-0.02012x2;y5=-50.75309+51.53527x-19.65947x2+3.66585x3-0.32886x4+0.01137x5;分别作出它们地曲线图,图中点划线为y1曲线,实线为y2曲线,虚线为y5曲线.’x’为给定地数据点.从图中可以看出并不是多项式次数越高越好,次数高了,曲线越能给定点处和实际吻合,但别地地方就很差了.因此,本例选用一次和两次地多项式拟合应该就可以了.实验报告四题目: Romberg 积分法摘要:对于实际地工程积分问题,很难应用Newton-Leibnitz 公式去求解.因此应用数值方法进行求解积分问题已经有着很广泛地应用,本文基于Romberg 积分法来解决一类积分问题.前言:(目地和意义)1. 理解和掌握Romberg 积分法地原理;2. 学会使用Romberg 积分法;3. 明确Romberg 积分法地收敛速度及应用时容易出现地问题. 数学原理:考虑积分⎰=ba dx x f f I )()(,欲求其近似值,通常有复化地梯形公式、Simpsion公式和Cotes 公式.但是给定一个精度,这些公式达到要求地速度很缓慢.如何提高收敛速度,自然是人们极为关心地课题.为此,记T 1,k 为将区间[a,b ]进行2k 等分地复化地梯形公式计算结果,记T 2,k 为将区间[a,b ]进行2k 等分地复化地Simpsion 公式计算结果,记T 3,k 为将区间[a,b ]进行2k 等分地复化地Cotes 公式计算结果.根据Richardson 外推加速方法,可以得到收敛速度较快地Romberg 积分法.其具体地计算公式为: 1. 准备初值,计算)]()([21,1b f a f ba T +-=2. 按梯形公式地递推关系,计算∑-=-+-+-+-+=1201,11,11))5.0(2(221k i k k k k i ab a f a b T T 3. 按Romberg 积分公式计算加速值1441,11,11,--=----+---m mk m m k m m m k m T T T m=2,…,k4. 精度控制.对给定地精度R ,若R T T m m <--1,11,则终止计算,并取1,m T 为所求结果;否则返回2重复计算,直至满足要求地精度为止. 程序设计:本实验采用Matlab 地M 文件编写.其中待积分地函数写成function 地方式,例如如下function yy=f(x,y); yy=x.^3;写成如上形式即可,下面给出主程序Romberg 积分法源程序%%% Romberg 积分法 clear%%%积分区间 b=3; a=1;%%%精度要求 R=1e-5;%%%应用梯形公式准备初值 T(1,1)=(b-a)*(f(b)+f(a))/2; T(1,2)=T(1,1)/2+(b-a)/2*f((b+a)/2); T(2,1)=(4*T(1,2)-T(1,1))/(4-1); j=2; m=2;%%%主程序体%%%while(abs(T(m,1)-T(m-1,1))>R);%%%精度控制 j=j+1; s=0;for p=1:2^(j-2);s=s+f(a+(2*p-1)*h/(2^(j-1))); endT(1,j)=T(1,j-1)/2+h*s/(2^(j-1)); %%%梯形公式应用 for m=2:j; k=(j-m+1);T(m,k)=((4^(m-1))*T(m-1,k+1)-T(m-1,k))/(4^(m-1)-1); end end%%%给出 Romberg 积分法地函数表 I=T(m,1)结果分析和讨论: 1. 求积分dx x10063.精确解I= 24999676.运行程序得Romberg 积分法地函数表为1.0e+007 *4.70101520000000 3.05022950000000 2.63753307500000 2.49996760000000 2.49996760000000 0 2.49996760000000 0 0由函数表知Romberg 积分给出地结果为2.4999676*10^7,与精确没有误差,精度很高.2. 求积分dx xx⎰10sin . 直接按前面方法进行积分,会发现系统报错,出现了0为除数地现象.出现这种情况地原因就是当x=0时,被积函数分母出现0,如果用一个适当地小数ε(最好不要小于程序给定地最小误差值,但不能小于机器地最大精度)来代替可以避免这个问题.本实验取R =ε,可得函数表为:0.92073548319659 0.93979327500190 0.94451351171417 0.94569085359489 0.94598501993743 0.94614587227034 0.94608692395160 0.94608330088846 0.94608307538495 0 0.94608299406368 0.94608305935092 0.94608306035138 0 0 0.94608306038722 0.94608306036726 0 0 0 0.94608306036718 0 0 0 0故该函数地积分为0.94608306036718,取8位有效数字.3. 求积分dx x ⎰12sin本题地解析解很难给出,但运用Romberg 积分可以很容易给出近似解,函数表为:0.42073549240395 0.33406972582924 0.31597536075922 0.31168023948094 0.31062036680949 0.31035626065456 0.30518113697100 0.30994390573588 0.31024853238818 0.31026707591900 0.31026822526959 0 0.31026142365354 0.31026884083167 0.31026831215439 0.31026830189296 0 0 0.31026895856465 0.31026830376269 0.31026830173008 0 0 0 0.31026830119484 0.31026830172211 0 0 0 0 0.31026830172262 0 0 0 0 0故该函数地积分为0.31026830172262,取8位有效数字.结论:Romberg 积分通常要求被积函数在积分区间上没有奇点.如有奇点,且奇点为第一间断点,那么采用例3地方法,还是能够求出来地,否则,必须采用其它地积分方法.当然,Romberg 积分地收敛速度还是比较快地.。
数值分析上机实践报告
数值分析上机实践报告班级:计算机1002姓名:陈斯琪学号:20102686课题三A . 实验题目:线性方程组的迭代法B . 实验要求(1) 应用迭代法求解线性方程组,并与直接法作比较;(2) 分别对不同精度要求,如5-4-3-10,10,10=ε,利用所需迭代次数体会该迭代法的收敛快慢;(3) 对方程组(2),(3)使用SOR 方法时,选取松弛因子=0.8,0.9,1,1.1,1.2等,试观察对算法收敛性的影响,并找出你所选用松弛因子的最佳值;(4) 编制出各种迭代法的程序并给出计算结果。
C . 目的和意义(1) 通过上机了解迭代法求解线性方程组的特点;掌握求解线性方程组的各类迭代法;(2) 体会上机计算时,终止准则‖X^(k+1)-X^k ‖∞<ε,对控制迭代精度的有效性; (3) 体会初始值和松弛因子的选择,对迭代收敛速度的影响 D . 实验方程组(1)线性方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡1-421534100368-24-3-81-012029137-2621-234179-11-1003524-31-23-6217758-6233-761-62911-31-512-301-231-2-2010563-5-6000121-3-2416084-0484⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡10987654321x x x x x x x x x x =⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡-2119381346323125精确解T x )2,1,1,3,0,2,1,0,1,1(*--=.(2) 对称正定线性方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡---=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡45152211236601924-3-360024-3-36014110-3-5211144-3-310-4221-8-13-4-1-612-53-8-1141-2312-1-204204-2004204-2487654321x x x x x x x x 精确解T *)2,0,1,1,2,0,1,1(--=x . (3)三对角线性方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡554141262135741-000000001-000000041-0000001-41-0000001-41-0000001-41-0000001-41-0000001-41-0000001-41-0000001-400000001-000000001-410987654321x x x x x x x x x x 精确解T x )1,1,0,3,2,1,0,3,1,2(*---=.E . 实验程序代码及截图(1) 应用Jacobi 迭代法求解方程组 代码如下:#include<iostream.h> #include<math.h>#define N 10 //十阶矩阵static double A[N][N]={4,2,-3,-1,2,1,0,0,0,0,8,6,-5,-3,6,5,0,1,0,0,4,2,-2,-1,3,2,-1,0,3,1,0,-2,1,5,-1,3,-1,1,9,4,-4,2,6,-1,6,7,-3,3,2,3,8,6,-8,5,7,17,2,6,-3,5,0,2,-1,3,-4,2,5,3,0,1,16,10,-11,-9,17,34,2,-1,2,2,4,6,2,-7,13,9,2,0,12,4,0,0,-1,8,-3,-24,-8,6,3,-1};//方程组左侧系数矩阵 static double B[N]={5,12,3,2,3,46,13,38,19,-21}; //右侧值。
数值分析上机实验
数值分析上机实验⽬录1 绪论 (1)2 实验题⽬(⼀) (2)2.1 题⽬要求 (2)2.2 NEWTON插值多项式 (3)2.3 数据分析 (4)2.3.1 NEWTON插值多项式数据分析 (4)2.3.2 NEWTON插值多项式数据分析 (6)2.4 问答题 (6)2.5 总结 (7)3 实验题⽬(⼆) (8)3.1 题⽬要求 (8)3.2 ⾼斯-塞德尔迭代法 (8)3.3 ⾼斯-塞德尔改进法—松弛法 (9)3.4 松弛法的程序设计与分析 (9)3.4.1 算法实现 (9)3.4.2 运算结果 (9)3.4.3 数据分析 (11)4 实验题⽬(三) (13)4.1 题⽬要求 (13)4.2 RUNGE-KUTTA 4阶算法 (13)4.3 RUNGE-KUTTA 4阶算法运算结果及数值分析 (14)总结 (16)附录A (17)1绪论数值分析是计算数学的⼀个主要部分,它主要研究各类数学问题的数值解法,以及分析所⽤数值解法在理论上的合理性。
实际⼯程中的数学问题⾮常复杂,所以往往需要借助计算机进⾏计算。
运⽤数值分析解决问题的过程:分析实际问题,构建数学模型,运⽤数值计算⽅法,进⾏程序设计,最后上机计算求出结果。
数值分析这门学科具有⾯向计算机、可靠的理论分析、好的计算复杂性、数值实验、对算法进⾏误差分析等特点。
本学期开设了数值分析课程,该课程讲授了数值分析绪论、⾮线性⽅程的求解、线性⽅程组的直接接法、线性⽅程组的迭代法、插值法、函数逼近与曲线拟合、数值积分和数值微分、常微分⽅程初值问题的数值解法等内容。
其为我们解决实际数学问题提供了理论基础,同时我们也发现课程中很多问题的求解必须借助计算机运算,⼈⼯计算量太⼤甚⾄⽆法操作。
所以学好数值分析的关键是要加强上机操作,即利⽤计算机程序语⾔实现数值分析的算法。
本报告就是基于此⽬的完成的。
本上机实验是通过⽤计算机来解答数值分析问题的过程,所⽤的计算⼯具是⽐较成熟的数学软件MATLAB。
数值分析上机实践报告
数值分析上机实践报告班级:计算机1002姓名:陈斯琪学号:20102686课题三A . 实验题目:线性方程组的迭代法B . 实验要求(1) 应用迭代法求解线性方程组,并与直接法作比较;(2) 分别对不同精度要求,如5-4-3-10,10,10=ε,利用所需迭代次数体会该迭代法的收敛快慢;(3) 对方程组(2),(3)使用SOR 方法时,选取松弛因子=0.8,0.9,1,1.1,1.2等,试观察对算法收敛性的影响,并找出你所选用松弛因子的最佳值; (4) 编制出各种迭代法的程序并给出计算结果。
C . 目的和意义(1) 通过上机了解迭代法求解线性方程组的特点;掌握求解线性方程组的各类迭代法;(2) 体会上机计算时,终止准则‖X^(k+1)-X^k ‖∞<ε,对控制迭代精度的有效性; (3) 体会初始值和松弛因子的选择,对迭代收敛速度的影响 D . 实验方程组(1)线性方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡1-421534100368-24-3-81-012029137-2621-234179-11-1003524-31-23-6217758-6233-761-62911-31-512-301-231-2-2010563-5-6000121-3-20416084-0484⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡10987654321x x x x x x x x x x =⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡-2119381346323125精确解Tx )2,1,1,3,0,2,1,0,1,1(*--=.(2) 对称正定线性方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡---=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡45152211236601924-3-360024-3-36014110-3-5211144-3-310-4221-8-13-4-1-612-53-8-1141-2312-1-204204-2004204-2487654321x x x x x x x x精确解T*)2,0,1,1,2,0,1,1(--=x .(3)三对角线性方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡554141262135741-000000001-0041-0000001-41-0000001-41-0000001-41-0000001-41-0000001-41-0000001-41-0000001-400000001-000000001-410987654321x x xx x x x x x x精确解Tx )1,1,0,3,2,1,0,3,1,2(*---=.E . 实验程序代码及截图(1) 应用Jacobi 迭代法求解方程组代码如下: #include<iostream.h> #include<math.h>#define N 10 //十阶矩阵 staticdoubleA[N][N]={4,2,-3,-1,2,1,0,0,0,0,8,6,-5,-3,6,5,0,1,0,0,4,2,-2,-1,3,2,-1,0,3,1,0,-2,1,5,-1,3,-1,1,9,4,-4,2,6,-1,6,7,-3,3,2,3,8,6,-8,5,7,17,2,6,-3,5,0,2,-1,3,-4,2,5,3,0,1,16,10,-11,-9,17,34,2,-1,2,2,4,6,2,-7,13,9,2,0,12,4,0,0,-1,8,-3,-24,-8,6,3,-1};//方程组左侧系数矩阵 static double B[N]={5,12,3,2,3,46,13,38,19,-21}; //右侧值static double Y[N]; //输出比较项static double Y[N];static double X[N]; //输出项static double G[N]; //X = BX' + G的G矩阵int i,j,k; //计数器double eps;int M=100;bool distance(){ //求两输出项的差的范数是否满足精度要求double temp=0;for (i=0;i<N;i++){temp=temp+fabs(X[i]-Y[i]);}if (temp>eps)return false;elsereturn true; //满足精度要求则结束程序}void main(){cout<<"最大迭代次数为100次"<<endl;cout<<"你希望的精度是多少?"<<endl;cout<<"eps=";cin>>eps;//形成迭代矩阵B,存放到A中for (i=0;i<N;i++){if (fabs(A[i][i])<eps){cout <<"打印失败"<<endl;return;}double T=A[i][i];for (j=0;j<N;j++){A[i][j]=-A[i][j]/T;}A[i][i] = 0;G[i]=B[i]/T;}int counter=0;while (counter<M){//迭代for (i=0;i<N;i++){double temp=0;for (j=0;j<N;j++){temp=temp+A[i][j]*Y[j];}X[i]=G[i]+temp;}if (distance()==true)break;else{//交换X,Y向量;for(i=0;i<N;i++){Y[i]=X[i];}}counter++;}//打印Xcout << "迭代次数为:"<<counter<<"次。
《数值分析》上机实验报告
数值分析上机实验报告《数值分析》上机实验报告1.用Newton 法求方程 X 7-X 4+14=0在(0.1,1.9)中的近似根(初始近似值取为区间端点,迭代6次或误差小于0.00001)。
1.1 理论依据:设函数在有限区间[a ,b]上二阶导数存在,且满足条件{}αϕ上的惟一解在区间平方收敛于方程所生的迭代序列迭代过程由则对任意初始近似值达到的一个中使是其中上不变号在区间],[0)(3,2,1,0,)(')()(],,[x |))(),((|,|,)(||)(|.4;0)(.3],[)(.20)()(.110......b a x f x k x f x f x x x Newton b a b f a f mir b a c x f ab c f x f b a x f b f x f k k k k k k ==-==∈≤-≠>+令)9.1()9.1(0)8(4233642)(0)16(71127)(0)9.1(,0)1.0(,1428)(3225333647>⋅''<-=-=''<-=-='<>+-=f f x x x x x f x x x x x f f f x x x f故以1.9为起点⎪⎩⎪⎨⎧='-=+9.1)()(01x x f x f x x k k k k 如此一次一次的迭代,逼近x 的真实根。
当前后两个的差<=ε时,就认为求出了近似的根。
本程序用Newton 法求代数方程(最高次数不大于10)在(a,b )区间的根。
1.2 C语言程序原代码:#include<stdio.h>#include<math.h>main(){double x2,f,f1;double x1=1.9; //取初值为1.9do{x2=x1;f=pow(x2,7)-28*pow(x2,4)+14;f1=7*pow(x2,6)-4*28*pow(x2,3);x1=x2-f/f1;}while(fabs(x1-x2)>=0.00001||x1<0.1); //限制循环次数printf("计算结果:x=%f\n",x1);}1.3 运行结果:1.4 MATLAB上机程序function y=Newton(f,df,x0,eps,M)d=0;for k=1:Mif feval(df,x0)==0d=2;breakelsex1=x0-feval(f,x0)/feval(df,x0);ende=abs(x1-x0);x0=x1;if e<=eps&&abs(feval(f,x1))<=epsd=1;breakendendif d==1y=x1;elseif d==0y='迭代M次失败';elsey= '奇异'endfunction y=df(x)y=7*x^6-28*4*x^3;Endfunction y=f(x)y=x^7-28*x^4+14;End>> x0=1.9;>> eps=0.00001;>> M=100;>> x=Newton('f','df',x0,eps,M);>> vpa(x,7)1.5 问题讨论: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*的迭代数列{x k},这就是Newton法的思想。
当x0接近x*时收敛很快,但是当x0选择不好时,可能会发散,因此初值的选取很重要。
另外,若将该迭代公式改进为其中r为要求的方程的根的重数,这就是改进的Newton法,当求解已知重数的方程的根时,在同种条件下其收敛速度要比Newton法快的多。
程序设计:本实验采用Matlab的M文件编写。
其中待求解的方程写成function的方式,如下function y=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;if f12(a)*f12(c)>0;a=c;elseb=c;endR=b-a;%求出误差k=k+1;endx=c%给出解Newton法及改进的Newton法源程序:clear%%%% 输入函数f=input('请输入需要求解函数>>','s')%%%求解f(x)的导数df=diff(f);%%%改进常数或重根数miu=2;%%%初始值x0x0=input('input initial value x0>>');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);breakendif k>max;%如果迭代次数大于给定值,认为迭代不收敛,重新输入初值ss=input('maybe result is error,choose a new x0,y/n?>>','s');if strcmp(ss,'y')x0=input('input initial value x0>>');k=0;elsebreakendendendk;%给出迭代次数x=x0;%给出解结果分析和讨论:1.用二分法计算方程在[1,2]内的根。
哈工大数值分析实验报告
哈工大数值分析实验报告标题:哈工大数值分析实验报告一、实验目的:本实验的目的是探究在数值分析中使用的各种数值方法,对于解决实际问题的有效性和可靠性进行评估。
二、实验内容:本实验主要包括以下几个方面的内容:1. 熟悉数值分析中常用的数值方法,如数值积分、数值微分、迭代法等;2. 在MATLAB等数学软件平台上,编写程序实现所学的数值方法;3. 使用所编写的程序,对给定的实际问题进行求解,并分析其结果的有效性和可靠性;4. 根据实际问题的特点,评估不同数值方法的适用性,并给出相应的结论和建议。
三、实验步骤:1. 阅读相关的理论知识,熟悉数值分析中常用的数值方法;2. 编写数值分析实验的程序代码,包括数值积分、数值微分和迭代法等;3. 使用编写的程序,对所给的实际问题进行求解,记录并分析结果;4. 根据实际问题的特点,评估所使用的数值方法的可靠性和有效性;5. 根据实验结果,撰写实验报告,包括实验目的、实验内容、实验步骤和实验结果的分析等。
四、实验结果:根据实际问题的不同,实验结果也会有所差异。
在实验报告中,可以详细叙述对所给实际问题的求解过程,并对结果进行分析和解释。
同时,还可以比较不同数值方法的结果,评估其优劣和适用性。
五、实验结论:根据实验结果的分析,可以得出结论,总结不同数值方法的优缺点,并对其在实际问题中的应用进行评价。
同时,还可以给出相应的建议,为以后的数值分析工作提供参考。
六、实验总结:通过本次实验,进一步加深了对数值分析中常用数值方法的理解和掌握。
通过实际问题的求解,对于这些数值方法的应用和效果有了更深入的认识。
同时,也提高了编程和科研报告撰写的能力,为以后的学习和工作打下了坚实的基础。
以上是关于哈工大数值分析实验报告的基本内容,具体实验细节和结果请根据实际情况进行补充。
哈工大-数值分析上机实验报告
Emax= 0.70770085900503,0 此时由 Emax 可以看出,不选主元的结果应该可以说是不正确了,这是由机器误差引 起的。 当 10 20 时,不选主元和选主元的计算结果如下 NaN NaN NaN Emax=NaN, 0 不选主元时,程序报错: Warning: Divide by zero. 。这是因为机器计算的最小精度为 10-15,所以此时的 10 20 就认为是 0,故出现了错误现象。而选主元时则没有这种现象, 而且由 Emax 可以看出选主元时的结果应该是精确解。
x3 x 1 0
x0=1; x0=0.45, x0=0.65;
( x 1) 2 (2 x 1) 0
当 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 法收敛与初值很有关系, 有的时候甚至可能不收敛。
实验报告
结果分析和讨论: 例 用最小二乘法处理下面的实验数据 . xi fi 3 2.01 4 2.98 5 3.50 6 5.02 7 5.47 8 6.02 9 7.05
Hale Waihona Puke 并作出 f ( x) 的近似分布图。 分别采用一次,二次和五次多项式来拟合数据得到相应的拟合多项式为: y1=-0.38643+0.82750x ; y2=-1.03024+1.06893x-0.02012x2; y5=-50.75309+51.53527x-19.65947x2+3.66585x3-0.32886x4+0.01137x5; 分别作出它们的曲线图,图中点划线为 y1 曲线,实线为 y2 曲线,虚线为 y5 曲线。’x’为 给定的数据点。从图中可以看出并不是多项式次数越高越好,次数高了,曲线越能给定点 处和实际吻合,但别的地方就很差了。因此,本例选用一次和两次的多项式拟合应该就可 以了。
《数值分析》上机实验报告
数值分析上机实验报告x k x k - f(X k) f (X k)《数值分析》上机实验报告1. 用Newt on法求方程X7-X4+14=0在(0.1,1.9)中的近似根(初始近似值取为区间端点,迭代6次或误差小于0.00001 )。
1.1理论依据:设函数在有限区间[a,b]上二阶导数存在,且满足条件1. f(x)f(b) 02. f(x)在区间[a, b]上不变号3f(x) = 0;4」f (c)〔f .(x) |,其中c是a,b中使mir(| f .(a), f .(b) |)达到的一个b -a则对任意初始近似值x0• [a,b],由Newton迭代过程込f(x k )X“ M(Xk) = Xk — T^,k = 0,1,2,3…f'(X k)所生的迭代序列 % [平方收敛于方程f(x)=0在区间[a,b]上的惟一解: 令7 4f(x)=x -28x 14, f (0.1) 0, f(1.9) ::0f (x) =7x6-112x3=7x3(x3-16) ::: 0f (x) =42x5-336x2=42x2(x3-8) :: 0f (1.9) f (1.9) 0故以1.9为起点x0 =1.9如此一次一次的迭代,逼近X的真实根。
当前后两个的差<=出寸,就认为求出了近似的根。
本程序用Newton法求代数方程(最高次数不大于10)在(a,b )区间的根//限制循环次数1.2 C 语言程序原代码:#i nclude<stdio.h> #in clude<math.h> mai n() {double x2,f,f1; double x1=1.9; // 取初值为 1.9do {x2=x1;f=pow(x2,7)-28*pow(x2,4)+14; f1=7*pow(x2,6)-4*28*pow(x2,3); x 仁 x2-f/f1;}while(fabs(x1-x2)>=0.00001||x1<0.1); printf("计算结果:x=%f\n",x1);}1.3运行结果:* D:\VC + +\EXERCIS E\Debu g\l1.4 MATLAB上机程序fun cti on y=Newt on( f,df,x0,eps,M)d=0;for k=1:Mif feval(df,x0)==0d=2; breakelsex1=x0-feval(f,x0)/feval(df,x0);ende=abs(x1-x0);x0=x1;if e<=eps&&abs(feval(f,x1))v=epsd=1; breakendendif d==1y=x1;elseif d==0y='迭代M次失败';elsey=奇异’endfun cti on y=df(x)y=7*x A6-28*4*x A3;Endfunction y=f(x) y=x A7-28*x A4+14;End>> x0=1.9;>> eps=0.00001;>> M=100;>> x=Newto n('f,'df,x0,eps,M);>> vpa(x,7)1.5问题讨论:1•使用此方法求方解,用误差来控制循环迭代次数,可以在误差允许的范围内得到比较理想的计算结果。
数值分析上机实验——数值积分
实验报告哈尔滨工程大学教务处制实验三 数值积分一.数值积分的基本思想1.复合梯形公式:Tn=++)()([2b f a f h2∑-=11)](n k xk f ;2.复合辛普森公式:Sn=6h[f(a)+f(b)+2∑-=11)](n k xk f +4∑-=+1)2/1(n k x f ];以上两种算法都是将a-b 之间分成多个小区间(n ),则h=(b-a)/n,x k =a+kh,x k+1/2=a+(k+1/2)h,利用梯形求积根据两公式便可。
3.龙贝格算法:在指定区间内将步长依次二分的过程中运用如下公式(1)Sn=34T2n-31Tn(2)Cn=1516S2n-151Sn(3)Rn=6364C2n-631Cn4T)(k m=144-m m T )1(1+-k m - 141-mT )(1k m -,k = 1,2,… 二.实验题目及实验目的(第4章计算实习题第1题)用不同数值方法计算积分xdx x ln 1⎰= -94。
(1)取不同的步长h 。
分别用复合梯形及复合辛普森求积计算积分,给出误差中关于h 的函数,并与积分精确值比较两个公式的精度,是否存在一个最小的h ,使得精度不能再被改善?(2)用龙贝格求积计算完成问题(1)。
(3)用自适应辛普森积分,使其精度达到104-。
三.实验手段:指操作环境和平台:win7系统下MATLAB R2009a程序语言:一种类似C 语言的程序语言,但比C 语言要宽松得多,非常方便。
四.程序①复合梯形求积程序function t=TiXing_quad(a,b,.h) format longx=a:h:b;y=sqrt(x).*log(x);y(1)=0;t=0;for k=1:(b-a)/h,t= t+y(k)+y(k+1);endt=t*h/2;②复合辛普森求积程序function s=Simpson_quad(a,b,h) format longx=a:h:b;y=sqrt(x).*log(x);z=sqrt(x+h/2).*log(x+h/2);y(1)=0;s=0;for k=1:(b-a)/h,s= s+y(k)+y(k+1)+4*z(k);ends=s*h./6;③龙贝格求积程序function [q,R]=Romberg(a,b,eps) h=b-a;R(1,1)=h*(0+sqrt(b).*log(b))/2; M=1;J=0;err=1;while err>epsJ=J+1;h=h/2;S=0;for p=1:Mx=a+h*(2*p-1);S=S+sqrt(x).*log(x);endR(J+1,1)=R(J,1)/2+h*S;M=2*M;for k=1:JR(J+1,k+1)=R(J+1,k)+(R(J+1,k)-R(J,k))/(4^k-1);enderr=abs(R(J+1,J)-R(J+1,J+1));endq=R(J+1,J+1);控制台输入代码:(1)>> a=0;>> b=1;>> h=0.1;>> t=TiXing_quad(a,b,h)>> s=Simpson_quad(a,b,h)>> h=0.01;>> t=TiXing_quad(a,b,h)>> s=Simpson_quad(a,b,h)>> h=0.001;>> t=TiXing_quad(a,b,h)>> s=Simpson_quad(a,b,h)(2)>> a=0;>> b=1;>> eps=10^-8;>> [quad,R]=Romberg(a,b,eps)(3)>> a=0;>> b=1;>> eps=10^-4;>> q=ZiShiYingSimpson('sqrt(x).*log(x)',a,b,eps) 五.实验结果比较与分析(1)h = 0.1时h = 0.01时h = 0.001时由结果(1)可知对于同一步长h,复合辛普森法求积分精度明显比复合梯形法求积的精度要高,且当步长取不同值时即h越小时,积分精度越高。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
实验报告一题目:非线性方程求解摘要:非线性方程的解析解通常很难给出,因此线性方程的数值解法就尤为重要。
本实验采用两种常见的求解方法二分法和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*的迭代数列{x k},这就是Newton法的思想。
当x0接近x*时收敛很快,但是当x0选择不好时,可能会发散,因此初值的选取很重要。
另外,若将该迭代公式改进为其中r为要求的方程的根的重数,这就是改进的Newton法,当求解已知重数的方程的根时,在同种条件下其收敛速度要比Newton法快的多。
程序设计:本实验采用Matlab的M文件编写。
其中待求解的方程写成function的方式,如下function y=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;if f12(a)*f12(c)>0;a=c;elseb=c;endR=b-a;%求出误差k=k+1;endx=c%给出解Newton法及改进的Newton法源程序:clear%%%% 输入函数f=input('请输入需要求解函数>>','s')%%%求解f(x)的导数df=diff(f);%%%改进常数或重根数miu=2;%%%初始值x0x0=input('input initial value x0>>');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);breakendif k>max;%如果迭代次数大于给定值,认为迭代不收敛,重新输入初值ss=input('maybe result is error,choose a new x0,y/n?>>','s');if strcmp(ss,'y')x0=input('input initial value x0>>');k=0;elsebreakendendendk;%给出迭代次数x=x0;%给出解结果分析和讨论:1.用二分法计算方程在[1,2]内的根。
(,下同)计算结果为x= 1. 23;f(x)= -3. 311e-007;k=18;由f(x)知结果满足要求,但迭代次数比较多,方法收敛速度比较慢。
2.用二分法计算方程在[1,1.5]内的根。
计算结果为x= 1. 80;f(x)= 2. 815e-006;k=17;由f(x)知结果满足要求,但迭代次数还是比较多。
3.用Newton法求解下列方程a)x0=0.5;计算结果为x= 0. 78;f(x)= 2. 313e-016;k=4;由f(x)知结果满足要求,而且又迭代次数只有4次看出收敛速度很快。
b)x0=1;c)x0=0.45, x0=0.65;当x0=0.45时,计算结果为x= 0. 83;f(x)= -8. 584e-014;k=4;由f(x)知结果满足要求,而且又迭代次数只有4次看出收敛速度很快,实际上该方程确实有真解x=0.5。
当x0=0.65时,计算结果为x= 0. 00;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. 86;f(x)=4. 127e-007;k=16;显然这个结果不是很好,而且也不是收敛至方程的2重根上。
当x0=0.85时,结果收敛为x= 1. 89;f(x)= 2. 737e-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消去法源程序:cleara=input('输入系数阵:>>\n')b=input('输入列阵b:>>\n')n=length(b);A=[a b]x=zeros(n,1);%%%函数主体for k=1:n-1;%%%是否进行主元选取if abs(A(k,k))<yipusilong;%事先给定的认为有必要选主元的小数yzhuyuan=1;else yzhuyuan=0;endif yzhuyuan;%%%%选主元t=A(k,k);for r=k+1:n;if abs(A(r,k))>abs(t)p=r;else p=k;endend%%%交换元素if p~=k;for q=k:n+1;s=A(k,q);A(k,q)=A(p,q);A(p,q)=s;endendend%%%判断系数矩阵是否奇异或病态非常严重if abs(A(k,k))< yipusilongdisp(‘矩阵奇异,解可能不正确’)end%%%%计算消元,得三角阵for r=k+1:n;m=A(r,k)/A(k,k);for q=k:n+1;A(r,q)=A(r,q)-A(k,q)*m;endendend%%%%求解xx(n)=A(n,n+1)/A(n,n);for k=n-1:-1:1;s=0;for r=k+1:n;s=s+A(k,r)*x(r);endt=(A(k,n+1)-s)x(k)=(A(k,n+1)-s)/A(k,k)end结果分析和讨论:例:求解方程。
其中为一小数,当时,分别采用列主元和不列主元的Gauss消去法求解,并比较结果。
记E max为求出的解代入方程后的最大误差,按要求,计算结果如下:当时,不选主元和选主元的计算结果如下,其中前一列为不选主元结果,后一列为选主元结果,下同。
0. 91 0. 512. 72 2. 632. 51 2. 21E max= 9. 624e-010,0此时,由于不是很小,机器误差就不是很大,由E max可以看出不选主元的计算结果精度还可以,因此此时可以考虑不选主元以减少计算量。
当时,不选主元和选主元的计算结果如下1. 77 0. 481. 072. 743. 31 2. 09E max= 2. 668e-005,0此时由E max可以看出不选主元的计算精度就不好了,误差开始增大。
当时,不选主元和选主元的计算结果如下1. 20 1. 001. 662. 003. 11 000E max= 0. 03,0此时由E max可以看出,不选主元的结果应该可以说是不正确了,这是由机器误差引起的。
当时,不选主元和选主元的计算结果如下NaN 1NaN 2NaN 3E max=NaN, 0不选主元时,程序报错:Warning: Divide by zero.。
这是因为机器计算的最小精度为10-15,所以此时的就认为是0,故出现了错误现象。
而选主元时则没有这种现象,而且由E max可以看出选主元时的结果应该是精确解。
结论:采用Gauss消去法时,如果在消元时对角线上的元素始终较大(假如大于10-5),那么本方法不需要进行列主元计算,计算结果一般就可以达到要求,否则必须进行列主元这一步,以减少机器误差带来的影响,使方法得出的结果正确。