数值分析大作业三四五六七完整版.docx

1、数值分析大作业三四五六七完整版数值分析大作业三四五六七HEN system office room HEN 16H-HENS2AHENS8Q8-HENH1688大 作 业 三I.给定初值Xo及容许误差 ,编制牛顿法解方程fx) =0的通用 程序.解：Mat lab程序如下：函数m文件：function Fu=fu(x)Fu=x3/3-x;end函数m文件：function Fu=dfu(x)Fu=x*2-1;end用Newton法求根的通用程序clear;x0=inputC请输入初值x0:)；ep=input (，请输入容许i吴崔：)；flag=l;while flag=lxl=x0-fu(x

2、O)/dfu(xO);if abs(xl-xO)epflag=0;endx0=xl;endfprintf (方程的个近似解为:%fn, xO);寻找最大6值的程序：cleareps=input (请输入搜索淸度:;ep=input (，诘输入容许谋差：);flag=l;k=0;x0=0;while flag=lsigma=k*eps;x0=sigma;k=k+l;m=0;flagl=l;while flagl=l & m=103xl=xO-fu(xO)/dfu(xO);if abs(xl-xO)=epflag=O;endendfprintf 赧人的 sigma 值为：%fn*, sigma);

3、2.求下列方程的非零根解：Matlab程序为:（1）主程序clearclc format long x0=765;N=100;errorlim=10(-5);x=xO-f (xO)/subs (df (), xO); n=l;while nerrorlim n=n+l;elsebreak;endxO=x;enddisp( 迭代次数：n=*, num2str(n)J)disp(所求 IE零根:正根 xl=, num2str(x),负根 x2=, num2str (-x)J)(2)子函数非线性函数ffunction y=f(x)y=log(513+*x)/*x)-x/(1400*;end(3)函数

4、非线性函数的阶导数dffunction y=df()syms xly=log(513+*xl)/*xl)-xl/(1400*;y=diff(y);end运行结果如下：迭代次数:n=5所求非零根：正根xU负根x2=大作业四试编写 MATLAB函数实现Newton插值，要求能输出插值多项式.对函数/*(%)= 在区间-5, 5 1 + 4Q上实现10次多项式插值.分析：(1)输出插值多项式。 (2)在区间-5, 5内均匀插入99个节点，计算 这些节点上函数f (x)的近似值，并在同一张图上画出原函数和插值多项式的图 形。(3)观察龙格现象，计算插值函数在各节点处的误差，并画出误差图。解：Matla

5、b程序代码如下：%此函数实现y=l/(l+4*x*2)的n次Newton插值，n由调用函数时指定%函数输出为插值结果的系数向量(行向量和插值多项式function t y=func5(n)xO=linspace(-5, 5,n十1);y0=l./(l.+4. *x0. 2);b=zeros (1, n+1);for i=l:n+ls=0;for j=l:it=l；for k=l:iif k=jt=(xO(j)-xO(k)*t;end;end;s=s+yO(j)/t;end;b(i)=s;end;t=linspace(0,0, n+1);for i=l:ns=l inspace (0, 0, n

6、+1);s (n+l-i:n+1) =b(i+1). *poly(x0(l:i);t=t+s;end;t (n+l)=t (n+l)+b(l);y=poly2sym(t);10次插值运行结果：b Y=func5(10)b 二Columns 1 through 4Columns 5 through 8Columns 9 through 11Y 二b为插值多项式系数向量，Y为插值多项式。插值近似值：xl=linspace (-5, 5, 101);x=xl(2:100);y=polyval(b, x)y 二Columns 1 through 12绘制原函数和拟合多项式的图形代码: plot (x,

7、 1. / (1+4. *x2)hold allplot (x, y, r)xlabel ( X) ylabel ( Y) title ( Runge现象) gtextC原函数) gtext 十次牛顿插值多项式) 詁差计数并绘制误差图：hold offey=l. /(1+4 *x. 2)-yey =Columns 1 through 12Columns 13 through 24Columns 25 through 36Columns 37 through 48Columns 49 through 600Columns 61 through 72Columns 73 through 84Col

8、umns 85 through 96plot (x, ey)xlabelf X)ylabel ( ey)title ( Runge现象误差图)输出结果为：大作业五解：Mat 1 ab程序为：x = -520, -280, -7& , 0, 78, 280, 520;y = 0, -30, -36, -35, 0, 35, 36, 30, 0;n = 13;喘求解Mfor i = 1:l:n-lh(i) = x(i+l)-x(i);endfor i = 2:l:n-la(i) = h(i-l)/(h(i-l)+h(i);b(i) = l-a(i);c (i) = 6*(y(i+l)-y (i)/

9、h(i)-(y (i)-y(i-l)/h(i-l)/(h(i-l)+h(i);enda(n) = h(n-l)/(h(l)+h(n-l);b(n) = h(l)/ (h (l)+h(n-l);c(n) = 6/(h(l)+h(n-l)*(y-y(l)/h-(y(n)-y(n-l)/h(n-l);A(l, 1) = 2;A(l,2) = b(2);A(l, n-1) = a(2);A(n-1, n-2) = a(n);A(n-1, n-1) = 2;A (n-1,1) = b(n);for i = 2:1:n-2A(iti) = 2;A(i,i+1) = b(i+l);A(i,i-1) = a(

10、i+l); endC = c(2:n);m = AC;M(l) = m(n-l);M(2:n) = m;xx = -520:520;for i = 1:51for j = 1:l:n-lif x(j) = xx(i) & xx(i) x(j-l)break;endendyy(i) = M(j+l)*(xx(i)-x(j)厂3/(6*h(j)-M(j) *(xx(i)-x(j+1) *3/(6*h(j) + (y(j+l)- M(j+1)*h(j) *2/6)*(xx(i)-x(j)/h(j)-(y (j)-M(j) *h(j) 2/6)*(xx(i)-x(j+l)/h(j); end;for

11、i = 52:101yy(i) = -yy(102-i);end;for i = 1:50xx(i) = -xx (i);endplot (xx, yy);hold on;for i = 1:l:n/2x(i) = -x(i);endplot (x, y, bd)；titleC机翼外形曲线);输出结果：运行文件，得到2. (1)编制求第一型3次样条插值函数的通用程序；(2)已知汽车门曲线型值点的数据如下：解：(1) Matlab编制求第一型3次样条插值函数的通用程序：function Sx = Threch(X, Y, dyO, dyn)%X为输入变量x的数值喘Y为函数值y的数值瓯dyO为左端

12、阶导数值%dyn为右端阶导数值瓯Sx为输出的函数农达式n = length(X)l;d = zeros (n+1, 1);h = zeros(1, n-1);fl = zeros (1, n-1);f2 = zeros(l, n-2);for i = l:n 喘求函数的阶筮瀚h(i) = X(i+1)-X(i);fl(i) = (Y(i+l)-Y(i)/h(i);endfor i = 2:n 嗚求函数的二阶签商f2(i) = (f 1 (i)-fl (i-1)/(X(i+1) -X(i-1);d(i) = 6*f2(i);endd(l) = 6*(fl(l)-dy0)/h(l);d(n+l)

13、= 6*(dyn-fl(n-l)/h(nl);%赋初值A = zeros (n+1, n+1);B = zeros (1, n-1);C = zeros(1, n-1);for i = l:n-lB(i) = h(i)/(h(i)+h(i+l);C(i) = l-B(i);endA(l,2) = 1；A(n+l,n) = 1;for i = l:n+lA(i,i) = 2;endfor i = 2:nA(i,i-1) = B(i-l);A(ifi+1) = C(i-l);endM = Ad;syms x;for i = l:nSx(i) = collect(Y(i) + (fl(i)-(M(i

14、)/3+M(i+l)/6)*h(i)*(x-X(i)+M(i)/2*(x- X(i)*2+(M(i+l)-M(i)/(6*h(i)*(x-X(i)*3);digits(4);Sx(i) = vpa(Sx(i);endfor i = l:ndispC S(x)=);fprintf ( %s (%d, %d)n , char(Sx(i),X(i), X(i+1);endS = zeros(1,n);for i = l:nx = X十；S(i) = Y (i) + (f 1 (i) - (M (i) /3+M (i+1) /6) *h (i) * (x-X (i) +M (i) /2* (x-X (

15、i) *2+ (M (i+1)- M(i)/(6*h(i) *(x-X(i)厂 3;enddisp(* S(i+);dispC i X(i+ S(i+ );for i = l:nfprintf C %d %. 4f %. 4fn , i, X(i)+, S(i);endend输出结果： X = 0 1 2 3 4 5 6 7 8 9 10; Y = ; Threch (X, Y,S(x) =*x*x2- *x3 +(0,1)S(x) =*x*x2- *x3 +(1,2)S(x)二*x*x2- *x3 +(2, 3)S(x) =2一 *x- *x3 +(3, 4)S(x) =*x*x2+ *x3

16、 -(4, 5)S(x)二2一 *x- *x3 +(5, 6)S(x) =*x*x2+ *x3 -(6, 7)S(x) =2一 *x- *x3 +(7, 8)S(x) =*x*x2+ *x3 -(& 9)S(x)二*x*x2+ *x3 -(9, 10)S(i+1X(i+S(i+12345678910 ans- *x3 - *x2 + *x + ， - *x3 - *x2 + *x + , - *x3 - *x2 + *x + , - *x3 + *x2 - *x + , *x3 - *x2 + *x - , - *x3 + *x2 - *x + , *x3 - *x2 + *x - , - *

17、x3 + *x2 -*x + , *x3 - *x2 + *x - , *x3 - *x2 + *x 一 大作业六1、炼钢厂出钢时所用的圣刚睡的钢包，在使用过程中111于钢液及炉渣对包衬耐火 材料的侵蚀，使其容积不断增大，经试验，钢包的容积与相应的使用次数的数据如 下：(使用次数X,容积y)x y23f/?* 1/x5679iOx y9io14161719g选用双曲线l/y = a对使 用最乘法数据进行拟合。解：Mat lab程序如下：function a=nihehanshu()x0=2 3 5 6 7 9 10 11 12 14 16 17 19 20;y0= 1;A=zeros(2, 2

18、);B=zeros(2,1);a=zeros(2,1);x=l /x0;y=l. /y0;A(l, 1)=14;A(l, 2)=sum(x);A (2,1)=A(1,2);A(2, 2)=sum(x. *2);B(l)=sum(y);B(2)=sum(x. *y);a=AB;y=l. /(a(l)+a(2)*1. /x0); subplot (1, 2, 2);plot (xO, yO-y, bd-); titleC拟合曲线误差); subplot (1,2, 1);plot (xO, yO, go);hold on;x=2:20;y=l. /(a(l)+a(2)*1. /x); plot (

19、x, y, r*-);legend (散点,，拟介曲线图 1 / y=a (1) +a (2)*l/x); title(最小二乘法拟介曲线);试验所求的系数为：nihehanshuans =则拟合曲线为 丄= 0.00897327 6262446 + 0.00084169 0019705 -y %拟合曲线图、散点图、误差图如下：2、下面给出的是乌鲁木齐最近1个月早晨7:00左右（新媼时间）的天气预报所得 到的温度，按照数据找出任意次曲线拟合方程和它的图像。用MATLAB编程对上述 数据进行最小二乘拟合。2008年10月-11月26日112345678910温度91011121314131211

20、9天数11121314151617181920温度101112131412111098盍21222324252627282930温度78911976531解：Matlab的程序如下：x=l:l:30;y=9,10,11, 12, 13,14,13,12, 11, 9, 10,11,12,13,14, 12,11,10, 9, 8, 7, & 9,11, 9, 7, 6, 5, 3, 1;al=polyfit (x, 5% 3) %三次多项式拟合a2= polyfit (x, y, 9) %九次多项式拟合%a3= polyfit (x, y, 15) %十五次多项式拟合%bl= polyval(

21、al, x)b2= polyval(a2, x)b3= polyval(a3, x)rl= sum(y-bl). 2)滋三次多项式i吴差平方和r2= sum(y-b2). *2) %九次次多项式误差平方和r3= sum(y-b3). 2)弔十五次多项式误差平方和plot (x, y/*) %用*画出 x,y 图像%hold onplot (x, bl, r ) %用红色线画出x, bl图像%hold onplot (x, b2, g ) %用绿色线画出x, b2图像%hold onplot (x, b3, b:o)弘用蓝色o线画出x, b3图像%试验结果为：al =Columns 1 thro

22、ugh 2Columns a2 =ColumnsColumnsColumnsColumnsColumns a3 =ColumnsColumnsColumnsColumnsColumnsColumnsColumnsColumns bl =ColumnsColumnsColumnsColumnsColumnsColumnsColumnsColumnsColumnsColumnsColumnsColumnsColumnsColumnsColumns b2 =ColumnsColumnsColumnsColumnsColumnsColumnsColumnsColumnsColumnsColumnsCo

23、lumnsColumnsthrough 4through 2through 4 through 6 through 8 through 10through 2through 4through 6through 8 through 10through 12through 14through 16through 2through 4through 6through 8 through 10through 12through 14through 16through 18through 20through 22through 24through 26through 28through 30throug

24、h 2through 4through 6through 8 through 10through 12through 14through 16through 18through 20through 223135791357911131513579111315171921232527291357911131517192123through 24Columns25Columns27Columns29through through through2830b3 =Columns 1 through 2Columns3through 4Columns5through 6Columns7through 8

25、Columns9through 10Columns11through12Columns13through14Columns15through16Columns17through18Columns19through20Columns21through22Columns23through24Columns25through26Columns27through28Columns29through30rl = r2 = r3 =其中的最后图像为:大作业七用Romberg求积法计算积分的近似值，要求误差不超过5x 1旷解：Matlab程序为：玄被积函数m文件：function Fx=fx(x)Fx=l/

26、(l+100*x*x);end%Romberg求积法计算积分的通用程序function Romberg0clear;a=input (请输入积分下限:a=); b=input (请输入积分I.限：b=); eps=input (，请输入允许持度：eps=);%=计算 Tn=%function Tn=T(n)Tn=0;h=(b-a)/n;x=zeros (1, n+1);for k=l:n+lx(k)=a+(k-l)*h;endfor j=l:nTn=Tn+h*(fx(x(j)+fx(x(j+1)/2;endend%=计算 Sn=%function Sn=S(n)Sn=4/3*T(2*n)-l/3*T(n);end%=计算 Cn=%function Cn=C(n)Cn=16/15*S(2*n)-1/15*S(n);end%=计算 Rn=%function Rn=R(n)Rn=64/63*C(2*n)-1/63*C(n);end賀=d十算满足允许精度的Rn，并打印输!,=%i=l；flag=l；while flag=lif abs(R(2i)-R(2(i-l)/255epsflag=0