北科大工程数值上机实验.docx

1、北科大工程数值上机实验工程数值计算上机实验报告姓 名 丁恒 班 级 机自 1404学 号 41440001 2017年10月1 数据的插值与回归1.1 种插值算法的实践1.待定系数法多项式插值程序:图像:2拉格朗日多项式插值x012y1.00002.71837.3891程序： % 拉格朗日线性插值函数的构造比较clear;x0=0;y0=1;x1=1;y1=2.7183;x2=2;y2=7.3891;x=0:0.01:2;for i=1:length(x) ye(i)=exp(x(i); L(1)=(x1-x(i)*(x2-x(i)/(x1-x0)*(x2-x0); L(2)=(x0-x(i)

2、*(x2-x(i)/(x0-x1)*(x2-x1); L(3)=(x0-x(i)*(x1-x(i)/(x0-x2)*(x1-x2); y(i)=y0* L(1)+y1* L(2)+y2* L(3);endplot(x,y,k-,x,ye,r-,x0,y0,ro,x1,y1,ro);结果：3、牛顿多项式程序：clear;xi=0;40;80;120;160;200;240;280;320;360;yi=50.3;51.4;32.7;9.22;6.22;14.7;23.6;32.5;41.4;50.3;n=length(xi)B=zeros(n,n);B(1,:)=yi;for i=1:n-1 f

3、or j=1:n-i B(i+1,j)=(B(i,j+1)-B(i,j)/(xi(i+j)-xi(j); endend x=1:1:360;for i=1:length(x) ye(i)=log(x(i); y1(i)=B(1,1)+B(2,1)*(x(i)-xi(1); y2(i)=B(1,1)+B(2,1)*(x(i)-xi(1)+B(3,1)*(x(i)-xi(1)*(x(i)-xi(2); y3(i)=B(1,1)+B(2,1)*(x(i)-xi(1)+B(3,1)*(x(i)-xi(1)*(x(i)-xi(2)+B(4,1)*(x(i)-xi(1)*(x(i)-xi(2)*(x(i)

4、-xi(3); y4(i)=B(1,1)+B(2,1)*(x(i)-xi(1)+B(3,1)*(x(i)-xi(1)*(x(i)-xi(2)+B(4,1)*(x(i)-xi(1)*(x(i)-xi(2)*(x(i)-xi(3)+B(5,1)*(x(i)-xi(1)*(x(i)-xi(2)*(x(i)-xi(3)*(x(i)-xi(4);endplot(x,ye,k-,x,y1,r-,x,y2,b-,x,y3,g-,x,y4,c-,xi,yi,ro)结果：n =104、分段线性插值程序：clear;xi=0;40;80;120;160;200;240;280;320;360;yi=50.3;51

5、.4;32.7;9.22;6.22;14.7;23.6;32.5;41.4;50.3;% 主程序n=length(xi);B=zeros(n,n);B(1,:)=yi;for i=1:n-1 for j=1:n-i B(i+1,j)=(B(i,j+1)-B(i,j)/(xi(i+j)-xi(j); endendb=B(2:end,1);%x=-5:0.1:5;x=2:0.1:10; for j=1:length(x) Y(1,j)=B(1,1); p(1,j)=1;end%牛顿插值for i=1:length(b) for j=1:length(x) p(i+1,j)=p(i,j)*(x(j)

6、-xi(i); Y(i+1,j)=Y(i,j)+b(i)*p(i+1,j); endendplot(xi,yi,ko-,x,Y(4,:),r-);结果：5、分段三次厄米特插值插值程序：clear;N=10;for i=1:N+1 xi(i)=-5+10*(i-1)/N; yi(i)=xi(i)4; ypi(i)=4*xi(i)3;end n=length(xi);B=zeros(n,n);B(1,:)=yi;for i=1:n-1 for j=1:n-i B(i+1,j)=(B(i,j+1)-B(i,j)/(xi(i+j)-xi(j); endendb=B(2:end,1);x=-5:0.1:

7、5;for j=1:length(x) Y(1,j)=B(1,1); p(1,j)=1;end for i=1:length(b) for j=1:length(x) p(i+1,j)=p(i,j)*(x(j)-xi(i); Y(i+1,j)=Y(i,j)+b(i)*p(i+1,j); endendsyms rfor i=1:n-1 i a=xi(i); b=xi(i+1); A(i,1)=(1+2*(r-a)/(b-a)*(r-b)/(a-b)2; A(i,2)=(1+2*(r-b)/(a-b)*(r-a)/(b-a)2; A(i,3)=(r-a)*(r-b)/(a-b)2; A(i,4)=

8、(r-b)*(r-a)/(b-a)2; H(i,1)=A(i,1)*yi(i)+A(i,2)*yi(i+1)+A(i,3)*ypi(i)+A(i,4)*ypi(i+1);endHr=vpa(collect(H),8);for i=1:length(Hr) Ce(i,:)=sym2poly(Hr(i);end for i=1:length(x) ye(i)=x(i)4; for j=1:n-1 if x(i)=xi(j)& x(i)=xi(j)& x(i) XX = 117.1271 KK = 117.32652.2 MATLAB内部数值积分函数的应用分别用梯形公式和辛普森公式求函数y=4/(1

9、+x2)在0,1区间内的积分1、quad函数、Trapz函数y=4./(1+x.2)程序:clear;x=0:0.01:1;y=4./(1+x.2);% 梯形公式z1=trapz(x,y);% 辛普森公式z2=quad(4./(1+x.2),0,1);z1=vpa(z1,8)z2=vpa(z2,8)结果：z1 = 3.141576 z2 = 3.14159272.3 数值微分算法的实践1、对求差商程序：不同步长对差商求导clear;h=zeros(17,1);h(1)=0.05;a0=1;d0=1/(3*a0(2/3);n=length(h);for i=2:n h(i)=h(i-1)/5;e

10、ndfor i=1:n a(i)=a0-h(i); b(i)=a0+h(i); c(i)=(b(i)(1/3)-a(i)(1/3)/(b(i)-a(i); e(i)=c(i)-d0; er(i)=abs(e(i)/d0; cp(i)=1/(a(i)(2/3)+b(i)(2/3)+(a(i)*b(i)(1/3); ep(i)=cp(i)-d0; erp(i)=abs(ep(i)/d0; d(i)=d0;endsubplot(2,1,1);semilogx(h,c,r-,h,d,ko-);subplot(2,1,2);semilogx(h,er,b-); figure;subplot(2,1,1)

11、;semilogx(h,cp,r-,h,d,ko-);subplot(2,1,2);semilogx(h,erp,b-);结果：2、对求差商程序：% 不同步长对差商求导clear;h=zeros(17,1);h(1)=0.05;a0=1;d0=3*a02;n=length(h);for i=2:n h(i)=h(i-1)/5;endfor i=1:n a(i)=a0-h(i); b(i)=a0+h(i); c(i)=(b(i)3-a(i)3)/(b(i)-a(i); e(i)=c(i)-d0; er(i)=abs(e(i)/d0; cp(i)=a(i)2+b(i)2+a(i)*b(i); ep

12、(i)=cp(i)-d0; erp(i)=abs(ep(i)/d0; d(i)=d0;endsubplot(2,1,1);semilogx(h,c,r-,h,d,ko-);subplot(2,1,2);semilogx(h,er,b-); figure;subplot(2,1,1);semilogx(h,cp,r-,h,d,ko-);subplot(2,1,2);semilogx(h,erp,b-);结果：有上述结果知：若步长h越小，截断误差越小，结果越精确（1）当，会由于截断误差而造成失真（2）当，则不会造成失真2.4 MATLAB内部数值微分函数的应用diff函数的应用：比较不同方法求得f

13、(x)的导数程序：%数值微分函数的调用比较x=-3:0.1:3;f=inline(sqrt(x.3+2*x.2-x+12)+(x+5).(1/6)+5*x+2);g=inline(3*x.2+4*x-1)./(2*sqrt(x.3+2*x.2-x+12)+1/6./(x+5).(5/6)+5);p=polyfit(x,f(x),5);dp=polyder(p);dpx=polyval(dp,x);dx=diff(f(x,3.1)/0.1;gx=g(x);plot(x,g(x),r-,x,dpx,bo,x,dx,g-)结果：3 方程与方程组的求解3.1 线性方程组的求解1、克莱姆法程序：% 克莱

14、姆法求解线性方程组A=2 2 3;4 7 7;-2 4 5;b=3;1;-7;x=inv(A)*b结果：x = 2 -2 12、数值迭代法程序：%Jaboci迭代A=2 2 3;4 7 7;-2 4 5;b=3;1;-7;n=length(b);N=200;H=zeros(N,n);eb=1e-3for i=1:n x(1,i)=0;endfor i=1:N i for j=1:n for k=1:n if j=k H(i,j)=H(i,j)+A(j,k)*x(i,k); end end x(i+1,j)=(b(j)-H(i,j)/A(j,j); end e(i)=abs(x(i+1,1)-x

15、(i,1); if e(i)eb break endendx结果：i = 20x = 0 0 0 1.5000 0.1429 -1.4000 3.4571 0.6857 -0.9143 2.1857 -0.9184 -0.5657 3.2669 -0.5404 0.2090 1.7269 -1.9329 0.3391 2.9243 -1.1831 0.8371 1.4274 -2.3653 0.7162 2.7910 -1.3890 1.0632 1.2942 -2.5152 0.8276 2.7738 -1.4243 1.1298 1.2295 -2.5720 0.8489 2.7986 -

16、1.4086 1.1494 1.1845 -2.6058 0.8464 2.8362 -1.3804 1.1584 1.1427 -2.6363 0.8388 2.8781 -1.3489 1.1661 1.0998 -2.6679 0.8304 2.9223 -1.3159 1.1742 1.0546 -2.7012 0.82172.9687 -1.2815 1.1829结论：数值迭代不收敛。3.2 非线性方程的求解1、roots函数程序：% 非线性方程的求根p=1 2 3;roots(p)结果：ans = -1.0000 + 1.4142i -1.0000 - 1.4142i2、fsolv

17、e函数程序：function F= fun( x ) global p rF=p*sin(x)-r;end% 非线性方程的求根clear;global p rp=3;r=2;x=fsolve(fun,0);结果:x =0.72973、fzero函数程序：function y= fun1( x )y=exp(-x)-sin(pi/2*x);end% 非线性方程的求根clear;x=fzero(fun1,0);结果：x =0.44363.3非线性方程的求解1、程序：% 非线性方程组求解的应用clear; x,y=solve(x2+x*y+y-3=0,x2-4*x+3=0)结果：x = 1 3 y

18、= 1 -3/24 微分方程的求解初始条件：程序：sgm=10;b=2.666667;r=28;x0=6;y0=5;z0=5;T=20;dt=0.005;n=T/dt;t=0:dt:20;x(1)=x0;y(1)=y0;z(1)=z0;for i=1:n x(i+1)=x(i)+(-sgm*x(i)+sgm*y(i)*dt; y(i+1)=y(i)+(r*x(i)-y(i)-x(i)*z(i)*dt; z(i+1)=z(i)+(-b*z(i)+x(i)*y(i)*dt;endplot(x,y,r-,x,z,b-)结果：、 初始条件： 结果:4.1 R-K算法及实践程序：%四阶R-K的应用xi=0:0.2:1;yi(1)=1;n=length(xi);h=(xi(n)-x