1、工程數(shù)值計算上機實驗報告姓 名 丁恒 班 級 機自 1404學 號 2017年10月目 錄1 數(shù)據(jù)的插值與回歸31.1 種插值算法的實踐31.2 MATLAB內部插值函數(shù)的應用31.3 回歸算法的實踐31.4 MATLAB內部回歸函數(shù)的應用32 數(shù)值積分與數(shù)值微分42.1 數(shù)值積分算法的實踐42.2 MATLAB內部數(shù)值積分函數(shù)的應用42.3 數(shù)值微分算法的實踐42.4 MATLAB內部數(shù)值微分函數(shù)的應用43 方程與方程組的求解53.1 線性方程組的求解53.2 非線性方程的求解54 微分方程的求解64.1 R-K算法及實踐64.2 基于MATLAB內部ode函數(shù)的微分方程求解64.3 基于S

2、imulink的微分方程求解61 數(shù)據(jù)的插值與回歸1.1 種插值算法的實踐1. 待定系數(shù)法多項式插值程序:圖像:2拉格朗日多項式插值x012y1.00002.71837.3891程序： % 拉格朗日線性插值函數(shù)的構造比較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)*(x2-x(i)/(x0-x1)*(x2-x1); L(3)=(x0-x(i)*(x1-x

3、(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 for j=1:n-i B(i+1,j)=(B(i,j+1)-B(i,j)/(xi(i+j)

4、-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)-xi(3); y4(i)=B(1,1)+B(2,1)*(x(i)-xi(1)+B(3,1

5、)*(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.4;32.7;9.22;6.22;14.7;23.6;32.5;41.4;50.3;%

6、 主程序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)-xi(i); Y(i+1,j)=Y(i,j)+b(i)*p(i+1,j); enden

7、dplot(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:5;for j=1:length(x) Y(1,j)=B(1,1); p(1,j)=1

8、;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)=(r-b)*(r-a)/(b-a)2; H(i,1)=A(i,1)*yi(i)+A(i

9、,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內部數(shù)值積分函數(shù)的應用分別用梯形公式和辛普森公式求函數(shù)y=4/(1+x2)在0,1區(qū)間內的積分1、 quad函數(shù)、Trapz函數(shù)y=4./(1+x.2)

10、程序: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. z2 = 3.2.3 數(shù)值微分算法的實踐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;endfor i=1:n a(i)=a0-h(i); b(i)=a0+h(i); c(i)=(b(i)(1/3)

11、-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);semilogx(h,cp,r-,h,d,ko-);subplot(2,1,2);semilogx(h,er

12、p,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(i)=cp(i)-d0; erp(i)=abs(ep(i)/d0; d(i)=d0;endsubplot(

13、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內部數(shù)值微分函數(shù)的應用diff函數(shù)的應用：比較不同方法求得f(x)的導數(shù)程序：%數(shù)值微分函數(shù)的調用比較x=-3:0.1:3;f=inline(sqrt(x.3+2*x.

14、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、 克萊姆法程序：% 克萊姆法求解線性方程組A=2 2 3;4 7 7;-2 4 5;b=3;1;-7;x=inv(A)*b結果：x

15、 = 2 -2 12、 數(shù)值迭代法程序：%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(i,1); if e(i)eb break endendx結果：i = 20x = 0 0 0 1.5

16、000 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 -1.4086 1.1494 1.1845 -2.6058 0.8464 2.8362 -1.3804 1

17、.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結論：數(shù)值迭代不收斂。3.2 非線性方程的求解1、 roots函數(shù)程序：% 非線性方程的求根p=1 2 3;roots(p)結果：ans = -1.0000 + 1.4142i -1.0000 - 1.4142i2、 fsolve函數(shù)程序：function F= fun( x ) global p rF=p*sin(x)-r;

18、end% 非線性方程的求根clear;global p rp=3;r=2;x=fsolve(fun,0);結果:x =0.72973、 fzero函數(shù)程序：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 = 1 -3/24 微分方程的求解初始條件：程序：sgm=10;b=2.;r=28;x0=6;y0

19、=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)-xi(1)/(n-1);for i

20、=1:n k1(i)=yi(i)-2*xi(i)/yi(i); k2(i)=yi(i)+h/2*k1(i)-(2*xi(i)+h)/(yi(i)+h/2*k1(1); k3(i)=yi(i)+h/2*k2(i)-(2*xi(i)+h)/(yi(i)+h/2*k2(1); k4(i)=yi(i)+h*k3(i)-(2*xi(i)+2*h)/(yi(i)+h*k1(1); yi(i+1)=yi(i)+h/6*(k1(i)+2*k2(i)+2*k3(i)+k4(i);endyi=vpa(yi,7)結果：yi = 1.0, 1., 1., 1., 1., 1., 1. 4.2 基于MATLAB內部ode函數(shù)的微分方程求解程序:function F= fun3_4_