大连理工矩阵上机作业

1、第一题i.全世界石油产量（百万桶每H）如下表所示,确定并画出经过这些点的g阶多项 式，并使用该多项式估计缈网年的石油产量,龙格现黜在这个例子中出现了吗？以你的 观点，插值多项或是描述这些敕据好的模型吗？请解释一年百万桶每a便io.年百万桶每日fxlU6；效仍曾1999检的J 99568.008191f5698(/320(H74-4S7199772.024200274 065四8侬枷刻面76777Lagrange插值函数function y=lagrange(x0,y0,x);n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.

2、0;for j=1:nif j=kp=p*(z-x0(j)/(x0(k)-x0(j);endends=p*y0(k)+s;endy(i)=s;endx0=1:10;y0=67.052,68.008,69.803,72.024,73.400,72.063,74.669,74.487,74.065,76 .777;lagrange(x0,y0,17)ans=-1.9516e+12x=1:0.1:10;x=x;plot(x0,y0, r );hold onplot(x,y, k);legend(原函数,拟合函数)拟合图像如下B 5 7 57.s. 一；L: 用二口10W199619961907199

3、B19992000200120022003隼拟合函数出现了龙格现象，运用多项式进行插值拟合时，效果并不好，高次多项 式会因为误差的不断累积，导致龙格现象的发生。第二题2.用最小二泉法拟合题目J中的f衿数据点.拟合曲统为f”直践；留施物理,以 范三次曲慢，并计舁它们的均方误差.使用所捋则的物合曲践总计9用阱的产量.在 均方误星意义下H哪个拟合最好？function fun =nihe(n)m=67.052*10A6,68.008*10A6,69.803*10A6,72.024*10A6,73.400*10A6,72.063*10A6,74.669*10A6,74.487*10A6,74.065*

4、10A6,76.777*10A6;w=1,2,3,4,5,6,7,8,9,10;d1=0;d2=0;d3=0;y1=polyfit(m,w,1);y2=polyfit(m,w,2);y3=polyfit(m,w,3);y2=poly2sym(s2);y3=poly2sym(s3);y4=poly2sym(s4);f1=subs(y1,17);f2=subs(y2,17);f3=subs(y3,17);for p=1:10;d1=d1+(subs(y1,w(p)-m(p)A2;d2=d2+(subs(y2,w(p)-m(p)A2;d3=d3+(subs(y3,w(p)-m(p)A2;endd1=

5、sqrt(d1);d2=sqrt(d2);d3=sqrt(d3);fun=f1 f2 f3;d2 d3 d4;return ;结果拟合形式预计201件产量均力误差直线83382272.733087388.906抛物线74410590.912601531.997三次曲线99041977.162396531.259三次函数的均方误差最小，拟合的最好3.用WeuTtori法求解方程I*十* + 4 3 = D的根.初值选择坳=-0.7.迭代7步并与真 值k=1相比技，井列出如下表格函数function f=fun(x)syms aa=x;f=a*a*a+a*a+a-3;梯度函数function df

6、=dfun(x)df=3*x*x+2*x+1;Newtorftfunction result=didainewton(x0)k=0;xk=x0;xi=1;e0=abs(x0-xi);ek=e0;m=zeros(7,1);n=zeros(7,1);p=zeros(7,1);result=zeros(7,3);while k=sigmal p=x;x=B*p+f;n=n+1;endeval( x );N=4时R (B) =2.58211,迭代法不收敛 N=8时R (B) = 6.04221,迭代法不收敛Gauss-Seidelfunction x=gaussseidel(N) syms MM=N;

7、A=zeros(M);b=ones(M,1);for i=1:Mfor j=1:MA(i,j)=1/(i+j-1);end endp=zeros(M,1);L=-tril(A,-1);D=diag(diag(A);U=-triu(A,1);G=(inv(D-L)*U;f=(inv(D-L)*b;x=G*p+f;n=1;detal=1e-4;while abs(x-p)=detalp=x;x=G*p+f;n=n+1;end eval( x );结果N=4X=2.8527 -14.7133 -3.5388 26.7350 T-6.4475N=8X=-3.920228.0256-9.4820-41.

8、3197-34.080628.8462 65.3427 T共钝梯度法function result=cg(N)syms MM=N;x0=zeros(M,1);A=zeros(M);b=ones(M,1);for i=1:Mfor j=1:MA(i,j)=1/(i+卜1)； end end k=0;r0=b-A*x0;p0=r0;sigmal=1e-4;rk=r0;pk=p0;while dot(r0,r0)=sigmalak=(dot(r0,r0)/(dot(pk,A*pk);ri=r0;xk=x0+ak*pk;rk=r0-ak*A*pk;bk=(dot(rk,rk)/(dot(ri,ri);

9、pk=rk+bk*p0;p0=pk;r0=rk;x0=xk;k=k+1;endresult=x0;return ;结果 N=4 X=-4.0000 60.0000 -180.0000 140.0000 TN=8X= 4.5773 -72.4211 199.8408 -9.9954-174.8581 -171.8511 -10.8652269.4734 T第五题.编程计算三次样条满足5(0) =3.S=15(1) = 2.其中边界条件S“=function f=san(x1,y) n=length(x1);h=zeros(1,n-1);for p=1:n-1h(p)=x1(p+1)-x1(p);

10、 endfor p=1:n-2r(p)=h(p+1)/(h(p+1)+h(p);u(p+1)=h(p)/(h(p+1)+h(p);end u(1)=1;r(n-1)=1;c=2*ones(1,n);A=diag(r,-1)+diag(u,1)+diag(c);g=zeros(n,1);g(1)=3*(y(2)-y(1)/h(1);g(n)=3*(y(n)-y(n-1)/h(n-1);for p=1;n-2;g(p+1)=3*(y(p+2)-y(p+1)/h(p+1)*u(p+1)+r(p)*(y(p+1)-y(p)/h(p);endm=Ag;m=m;f=m;returnx1=0:4;y=1,3

11、,3,4,2;san(x1,y) ans =2.5000 1.0000 -0.5000 1.0000 -3.5000利用教科书p179的公式(5-35)function fun=sanci(m,h,xk,y)syms xn=length(m);l=length(xk); for i=1:n-1fun=(h(i)+2*(x-xk(i)*y(i)/h(i)+(x-xk(i)*m(i)*(x-xk(i+1)A2/h(i)A2+(h(i)-2*(x-xk(i+1)*y(i+1)/h(i)+(x-xk(i+1)*m(i+1)*(x-xk(i)A2/h(i)A2endend0=x1(9*x)/2 + 1)

12、*(x - 1)A2 - xA2*(5*x - 8)1=x2(7*x - 4)*(x - 2)A2 - (13*x)/2 - 16)*(x - 122=x3(11*x)/2 - 8)*(x - 3)A2 - (7*x - 25)*(x - 223=x4(9*x - 23)*(x - 4)A2 - (15*x)/2 - 32)*(x - 32第六题.取不同的初值用强截法求方程/+ 2产+ 13-100 = Q的实报.列表或者画图说 明欣鼓速度一函数function f=fun(x)syms a;a=x;f=a*a*a+2*a*a+10*a-100;end梯度函数function grad=gfun(x)syms aa=x;grad=3*a*a+4*a+10;迭代函数f unction result=xianjiefa(x0)sigmal=1e-6;maxi=500;m=zeros(500,1);d0=fun(x0)/gfun(x0);x1=x0-d0;k=1;xk=x1;m(1)=x1;while kmaxidk=fun