《数值计算技巧Matlab实题训练》由会员分享,可在线阅读,更多相关《数值计算技巧Matlab实题训练(15页珍藏版)》请在金锄头文库上搜索。
1、数值计算方法训练实习报告题 目:6A组院 系:上海电力学院数理学院专业年级:信息与计算科学专业2009级学生姓名: XX远 学号: 200924262011年7月8日第1题:含炭量与时间的关系在某冶炼过程中,钢的含炭量y与时间t的统计数据如下t 0510152025303540455055y 01.22.12.83.43.84.14.34.54.54.04.676647571824(1) 画出原始数据分布趋势图;(2) 用最小二乘法求钢的含炭量y与时间t的拟合曲线)=at + bt2 + ct3 ;(3) 打印出拟合曲线;(4) 另外选用y = atb进行拟合,比较二种拟合的效果。解:分析:使
2、用到曲线拟合的最小二乘法,对于拟合函数,尽量转化为可以方便 提炼出基函数的方程。在明确基函数的基础上,通过计算,得到各个系数,得到 法方程组(1),程序:function yuan(y)t=0:5:55;plot(t,y,*)legend(原始数据分布趋势图)运行结果:yuan(0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.5854.543.534.02 4.64)*米原始数据分布趋势图2.5210.5原始数据分布趋势(2),彳so:取基函数为:中0 = t 中1 = t2 中2 - t3:由基函数和y求法方程组的系数:(%,%)=肉 (土)叩。(?i
3、=1(阡号=幻1(七)气(七)i=1的 & ) = & (x )叩(X )222i 2ii=1(f & ) = y 叩(x )0i 0 ii=1:由这些系数,的。,气)=&。(七)气(七)i=1的产2)=& 1(七)叩23)i=1(%,%) = & O(P%(七)i=1(f & ) = y 叩 3) 1i 1 ii=1确定法方程组:(f &2) = y .%(七)i=1(甲,甲0)(平,甲)01(甲,甲2)(f 平。(f, %)(f, /)2(也,甲)(平,甲)1 1(也,甲2)(甲2,甲0)(平,甲)21(甲2,甲2):解这个法方程组:A - X = Bb,得到拟合函数:y = at + b
4、t2 + ct3 c程序:function a,b,c=xian(y0) t0=0:5:55;k1=t0;k2=t0.*t0;k3=t0.*t0.*t0;A=sum(k1.*k1)sum(k2.*k1)sum(k3.*k1);sum(k1.*k2)sum(k2.*k2)sum(k3.*k2);sum(k1.*k3) sum(k2.*k3) sum(k3.*k3);B=sum(k1.*y0);sum(k2.*y0);sum(k3.*y0);x=pinv(A)*B;a=x(1,1);b=x(2,1);c=x(3,1);t=0:55;y=a.*t+b.*t.”2+c.*t.”3;plot(t,y,,
5、-,) hold on plot(t0,y0,,*,)legend(y=a*t+b*t2+c*t3 拟合效果,真实值) 运行结果:a,b,c=xian(0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64) a =0.2657 b =-0.0053 c =3.5168e-005(3)拟合的图形,即上一题显示的图像图2拟合函数J = at + bt2 + ct3效果(4),用于这种非线性模型的拟合:把其化作线性:J=atb一两边同时取以e为底的对数一ln( j) = ln a + b - ln t:重复上面第二题的步骤进行,其中需要强调
6、的是(0,0)的点需要另外输入,因为ln(0)不存在,为了在同图出现,故对第二条拟合函数,取m = a n = b 程序:function m,n,a,b,c=fei(y2,y0)%y2=y0 除了 0 以外的数y1=log(y2);t1=5:5:55;n=length(t1);k1=ones(1,n);k2=log(t1);A=sum(k1.*k1) sum(k2.*k1);sum(k1.*k2) sum(k2.*k2);B=sum(k1.*y1);sum(k2.*y1);x=pinv(A)*B;m=exp(x(1,1);n=x(2,1);t=0:55;y=m*t.n;plot(t,y,,-
7、,)hold ona,b,c=xian(y0)plot(t,y,,-,)hold onplot(t1,y2,,*,,0,0,,*,)legend(,y=m*t.”n 拟合效果,y=a*t+b*t2+c*t3 拟合效果,真实值,)得到的拟合图像:1102030405060654320木0y=m*t. n拟合效果y=a*t+b*t 2+c*t 3 拟合效果真实值图3两种拟合函数拟合效果对比结论:在实际生活当中,不免需要对一组数据进行拟合,通过采用最佳的拟 合,找到一个近似的函数来研究数据的共性。通过这一道题目,发现不同的函数, 拟合效果差别也是蛮大的。第2题:特征值与特征向量用幕法求下列矩阵的主特
8、征值与相应的特征向量2 46 -3- 4 3一(1) A =3 9 15(2) A =-4 6 3124 16 36331解:利用幕法求矩阵A的主特征值与相应的特征向量,首先要给一个初始向量:T1 :定义一个和A行数一致的1列全一矩阵,即v。= 1M1nxl :v0= A.v0,为了方便计算,减少计算量,需要求出V1中按模最大的那个分量的值P,同时得到向量u =匕,由此可知u = v11 P001:重复第二步,得到V2 P2:计算P -P|,看其是否大于给定的误差,如果大于,则令V = V,P = P2 11212并从第三步开始重复;如果小于,则P为所求的按模最大的特征值,即主特征值,v2为其
9、对应的特征向量流程图:表1幂法求主特征值的流程图程序:function v2,p2=maxtr(A,err)n=size(A,2);v0=ones(n,1);v1=A*v0;p1=max(v1);u1=v1/p1;v2=A*u1;p2=max(v2);i=1;while i0if abs(p2-p1) f 2 (0 x 5)J (0) = 2取步长h=0.25,并画出0 x 5数值解的曲线。解:对于四阶龙格-库塔法,我们在求第k+1个点的值时,必须知道第k个点的 值,在通过中间的一些变量,使结果更加精确。:由题意可以确定初始值,即丁。丁 2*1= h - f (Xn,*)(2):分别计算出:v
10、 k2 = h - f (x + h2,七k = h - f (x + /, y3 n 2 n(3) : y = y + + (k + 2k + 2k + k ) x = x + h,由 n = 0开始,直至x = 5 n+1n 6 1234n+1n程序:function x,y=rk(x0,y0,a,b,h)%x0,y0 为方程的初值 x(1)=x0;y(1)=y0;n=(b-a)/h;for j=1:nk1二h*f(x(j),y(j);k2=h*f(x(j)+h/2,y(j)+k1/2);k3=h*f(x(j)+h/2,y(j)+k2/2);k4=h*f(x(j)+h,y(j)+k3);y
11、(j+1)=y(j)+1/6*(k1+2*k2+2*k3+k4);x(j+1)=x(j)+h;endplot(x,y,r:)hold onxx=0:0.2:5;yy=2./(xx.”2+1);plot(xx,yy,g-)legend(h=0.25的数值解,解析解)function z=f(x,y)z=-x*y“2;运行结果:x,y=rk(0,2,0,5,0.25)x =Columns 1 through 600.25000.50000.75001.00001.2500Columns 7through 121.50001.75002.00002.25002.50002.7500Columns 13through 183.00003.25003.50003.75004.00004.2500Columns 19through 214.50004.75005.0000y =Columns 1through 62.00001.88231.599