西南交大数值分析第二次大作业(可以运行)_西南交大第2次作业
西南交大数值分析第二次大作业(可以运行)由刀豆文库小编整理，希望给你工作、学习、生活带来方便，猜你可能喜欢“西南交大第2次作业”。
数值分析第二次大作业（1）用Lagrange插值法 程序：
function f=Lang(x,y,x0)symst;f=0;n=length(x);for(i=1:n)
l=y(i);for(j=1:i-1)
l=l*(t-x(j))/(x(i)-x(j));end;for(j=i+1:n)
l=l*(t-x(j))/(x(i)-x(j));end
f=f+l;simplify(f);
if(i==n)if(nargin==3)
f=subs(f,'t',x0);else
f=collect(f);
f=vpa(f,6);end end end
x=[1,2,3,-4,5];y=[2,48,272,1182,2262];t=[-1];disp('插值结果')yt=Lang(x,y,t)disp('插值多项式')yt=Lang(x,y)ezplot(yt,[-1,5]);运行结果：
插值结果： Yt= 12.0000 插值多项式：
yt =4.0*t^43.0*t + 2.0
（2）构造arctan x在[1,6]基于等距节点的10次插值多项式 程序：
function f=New(x,y,x0)symst;if(length(x)==length(y))
n=length(x);c(1:n)=0.0;else disp('xºÍyάÊý²»µÈ£¡');return;end f=y(1);y1=0;xx=linspace(x(1),x(n),(x(2)-x(1)));for(i=1:n-1)for(j=1:n-i)y1(j)=y(j+1)-y(j);end
c(i)=y1(1);
l=t;for(k=1:i-1)
l=l*(t-k);end;
f=f+c(i)*l/factorial(i);simplify(f);
y=y1;
if(i==n-1)if(nargin==3)
f=subs(f,'t',(x0-x(1))/(x(2)-x(1)));else
f=collect(f);
f=vpa(f,6);end end end
>>x=[1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6];y=[atan(1),atan(1.5),atan(2),atan(2.5),atan(3),atan(3.5),atan(4),atan(4.5),atan(5),atan(5.5),atan(6)];disp('插值多项式')yt=New(x,y)ezplot(yt,[1,6]);hold on ezplot('atan(t)',[1,6])grid on
运行结果： 插值多项式
yt = 1.34684*10^(-10)*t^100.00000104758*t^70.00176296*t^4 + 0.0125826*t^3-0.0640379*t^2 + 0.250468*t + 0.785398（1）用MATLAB自带spline函数用于进行三次样条插值 程序：
>>x=[-5,-4,-3,-2,-1,0,1,2,3,4,5];y=[0.03846,0.05882,0.10000,0.20000,0.50000,1.00000,0.50000,0.20000,0.10000,0.05882,0.03846];xi=linspace(-5,5)yi=spline(x,y,xi);plot(x,y,'rp',xi,yi);hold on;syms x fx=1/(1+x^2);ezplot(fx);grid on 运行结果：
由图可知，三次样条插值多项式图像与原函数图像基本一致。（2）取第一类边界条件，用三弯矩法编写MATLAB程序 程序
functionyi=cubic_spline(x,y,ydot,xi)n=length(x);ny=length(y);h=zeros(1,n);lambda=ones(1,n);mu=ones(1,n);M=zeros(n,1);d=zeros(n,1);for k=2:n h(k)=x(k)-x(k-1);end
for k=2:n-1
lambda(k)=h(k+1)/(h(k)+h(k+1));mu(k)=1-lambda(k);d(k)=6/(h(k)+h(k+1))...*((y(k+1)-y(k))/h(k+1)-(y(k)-y(k-1))/h(k));end
d(1)=6/h(2)*((y(2)-y(1))/h(2)-ydot(1));d(n)=6/h(n)*(ydot(2)-(y(n)-y(n-1))/h(n));A=diag(2*ones(1,n));for i=1:n-1 A(i,i+1)=lambda(i);A(i+1,i)=mu(i+1);end M=Ad;for k=2:n if x(k-1)
+M(k)/6/h(k)*(xi-x(k-1))^3...+1/h(k)*(y(k)-M(k)*h(k)^2/6)*(xi-x(k-1))...+1/h(k)*(y(k-1)-M(k-1)*h(k)^2/6)*(x(k)-xi);return;end end
>>a =-5;b = 5;n = 10;h =(b-a)/n;x=a:h:b;y=1./(1+x.^2);ydot=[-5/338,5/338];
xx=a:0.01:b;yy=1./(1+xx.^2);m = length(xx);z = zeros(1,m);for i=1:m z(i)=cubic_spline(x,y,ydot,xx(i));end
plot(x,y,'o',xx,yy,'k:',xx,z,'k-');运行结果：本题直接利用MATLAB自带的cftool曲线拟合工具箱完成。拟合结果：
（1）用多项式拟合：
Linear model Poly1:
f(x)= p1*x + p2 Coefficients(with 95% confidence bounds):
p1 =
-0.8685(-0.8815,-0.8556)
p2 =
30.61(30.48, 30.73)
Goodne of fit:
SSE: 0.1733
R-square: 0.9993
Adjusted R-square: 0.9993
RMSE: 0.1112（2）最小二乘法插值：
Linear interpolant:
f(x)= piecewise polynomial computed from p Coefficients:
p = coefficient structure
Goodne of fit:
SSE: 0
R-square: 1
Adjusted R-square: NaN
RMSE: NaN R-square(确定系数)SSR：Sum of squares of the regreion，即预测数据与原始数据均值之差的平方和，公式为：
SST：Total sum of squares，即原始数据和均值之差的平方和，公式为：
“确定系数”是定义为SSR和SST的比值，故
其实“确定系数”是通过数据的变化来表征一个拟合的好坏。由上面的表达式可以知道“确定系数”的正常取值范围为[0 1]，越接近1，表明方程的变量对y的解释能力越强，这个模型对数据拟合的也较好
4.（1）16点复合梯形公式： 程序：
function I =T_quad(x,y)n=length(x);m=length(y);h=(x(n)-x(1))/(n-1);a=[1 2*ones(1,n-2)1];I=h/2*sum(a.*y);End 运行结果：
(1)x=0:0.0625:1 y=exp(x)I=T_quad(x,y)x =
Columns 1 through 7
0
0.0625
0.1250
0.1875
0.2500
0.3125
0.3750
Columns 8 through 14
0.4375
0.5000
0.5625
0.6250
0.6875
0.7500
0.8125
Columns 15 through 17
0.8750
0.9375
1.0000 y =
Columns 1 through 7
1.0000
1.0645
1.1331
1.2062
1.2840
1.3668
1.4550
Columns 8 through 14
1.5488
1.6487
1.7551
1.8682
1.9887
2.1170
2.2535
Columns 15 through 17
2.3989
2.5536
2.7183 I =
1.7188
(2)x=1:0.0625:2 y=sin(x)/x I=T_quad(x,y)x =
Columns 1 through 7 1.0000
1.0625
1.1250
1.1875
1.2500
1.3125
1.3750
Columns 8 through 14
1.4375
1.5000
1.5625
1.6250
Columns 15 through 17
1.8750
1.9375
2.0000 y = 0.6116 I = 0.6116（2）8点的复合simpson公式： 程序：
function I=S_quad(x,y)n=length(x);m=length(y)N=(n-1)/2;h=(x(n)-x(1))/N;a=zeros(1,n);for k=1:N a(2*k-1)=a(2*k-1)+1;a(2*k)=a(2*k)+4;a(2*k+1)=a(2*k+1)+1 end I= h/6*sum(a.*y)end 运行结果：
（1）x=0:0.0625:1 y=exp(x)I=S_quad(x,y)I =
1.7183
（2）x=1:0.0625:2 y=sin(x)/x I=S_quad(x,y)I =
1.6875
1.7500 1.8125
0.6116
（3）三点的Gau-Legendre积分公式
程序：
function I=G_quad(fun,a,b,N)h=(b-a)/N;I=0;for k=1:N
t=[-sqrt(3/5)0 sqrt(3/5)];A=[5/9 8/9 5/9];
F=feval(fun,h/2*t+a+(k-1/2)*h);
I=I+sum(A.*F)end I=h/2*I;end
运行结果：
(1)fun=inline('exp(x)')I=G_quad(fun,0,1,4)I =
2.2722 I =
5.1898 I =
8.9360 I =
13.7463 I =
1.7183(2)fun=inline('sin(x)/x')I=G_quad(fun,1,2,4)I =
1.5954 I =
3.0143 I =
4.2363 I =
5.2479 I =
0.6560
5、给定初值问题y’=-1000(y-x^2)+2x，与y（0）=0，请分别用Euler和预测-矫正Euler算法按步长h=0.1，0.01，0.001，0.0001计算其数值，分析其中遇到的现象及问题
解：首先根据可计算出原函数为：y=Ce(-1000x)+x2，其中C为常量，为简便计算取C=1，计算初值x0=0，y0=0。（1）Euler算法 程序：
function[x,y]=Euler_f(ydot_fun,x0,y0,h,N)x=zeros(1,N+1);x(1)=x0;y(1)=y0;for n=1:N x(n+1)=x(n)+h;y(n+1)=y(n)+h*feval(ydot_fun,x(n),y(n));end
format short e ydot_fun=inline('-1000*(y-x*x)+2*x','x','y')[x,y]=Euler_f(ydot_fun,0,0,0.1,10)xx=0:0.1:1 yy=exp(-1000*xx)+xx.*xx d=abs(y-yy)plot(xx,d)title('Eular方法误差趋势图 h=0.1')运行结果： h=0.1 x =
Columns 1 through 6
0 1.0000e-001 2.0000e-001 3.0000e-001 4.0000e-001 5.0000e-001
Columns 7 through 11
6.0000e-001 7.0000e-001 8.0000e-001 9.0000e-001 1.0000e+000 y =
Columns 1 through 6
0
0 1.0200e+000-9.6940e+001 9.6061e+003-9.5099e+005
Columns 7 through 11
9.4148e+007-9.3207e+009 9.2274e+011-9.1352e+013 9.0438e+015 yy =
Columns 1 through 6
1.0000e+000 1.0000e-002 4.0000e-002 9.0000e-002 1.6000e-001 2.5000e-001
Columns 7 through 11
3.6000e-001 4.9000e-001 6.4000e-001 8.1000e-001 1.0000e+000
d =
Columns 1 through 6
1.0000e+000 1.0000e-002 9.8000e-001 9.7030e+001 9.6060e+003 9.5099e+005
Columns 7 through 11
9.4148e+007 9.3207e+009 9.2274e+011 9.1352e+013 9.0438e+015
h=0.01 x =
Columns 1 through 6 0 1.0000e-002 2.0000e-002 3.0000e-002 4.0000e-002 5.0000e-002
Columns 7 through 11
6.0000e-002 7.0000e-002 8.0000e-002 9.0000e-002 1.0000e-001 y =
Columns 1 through 6
0
0 1.2000e-003-6.4000e-003 6.7200e-002-5.8800e-001
Columns 7 through 11
5.3180e+000-4.7825e+001 4.3047e+002-3.8742e+003 3.4868e+004 xx =
Columns 1 through 6 0 1.0000e-002 2.0000e-002 3.0000e-002 4.0000e-002 5.0000e-002
Columns 7 through 11
6.0000e-002 7.0000e-002 8.0000e-002 9.0000e-002 1.0000e-001 yy =
Columns 1 through 6
1.0000e+000 1.4540e-004 4.0000e-004 9.0000e-004 1.6000e-003 2.5000e-003
Columns 7 through 11
3.6000e-003 4.9000e-003 6.4000e-003 8.1000e-003 1.0000e-002
d =
Columns 1 through 6
1.0000e+000 1.4540e-004 8.0000e-004 7.3000e-003 6.5600e-002 5.9050e-001
Columns 7 through 11
5.3144e+000 4.7830e+001 4.3047e+002 3.8742e+003 3.4868e+004
4.9580e+000 4.8315e+001 4.2983e+002 3.8750e+003 3.4867e+004
h=0.001 x =
Columns 1 through 6 0 1.0000e-003 2.0000e-003 3.0000e-003 4.0000e-003 5.0000e-003
Columns 7 through 11
6.0000e-003 7.0000e-003 8.0000e-003 9.0000e-003 1.0000e-002 y = Columns 1 through 6
0
0 3.0000e-006 8.0000e-006 1.5000e-005 2.4000e-005
Columns 7 through 11
3.5000e-005 4.8000e-005 6.3000e-005 8.0000e-005 9.9000e-005 xx =
Columns 1 through 6 0 1.0000e-003 2.0000e-003 3.0000e-003 4.0000e-003 5.0000e-003
Columns 7 through 11
6.0000e-003 7.0000e-003 8.0000e-003 9.0000e-003 1.0000e-002
yy =
Columns 1 through 6
1.0000e+000 3.6788e-001 1.3534e-001 4.9796e-002 1.8332e-002 6.7629e-003
Columns 7 through 11
2.5148e-003 9.6088e-004 3.9946e-004 2.0441e-004 1.4540e-004 d =
Columns 1 through 6
1.0000e+000 3.6788e-001 1.3534e-001 4.9788e-002 1.8317e-002 6.7389e-003
Columns 7 through 11
2.4798e-003 9.1288e-004 3.3646e-004 1.2441e-004 4.6400e-005
h=0.0001 x =
Columns 1 through 6 0 1.0000e-003 2.0000e-003 3.0000e-003 4.0000e-003 5.0000e-003
Columns 7 through 11
6.0000e-003 7.0000e-003 8.0000e-003 9.0000e-003 1.0000e-002 y =
Columns 1 through 6
0
0 3.0000e-006 8.0000e-006 1.5000e-005 2.4000e-005
Columns 7 through 11
3.5000e-005 4.8000e-005 6.3000e-005 8.0000e-005 9.9000e-005 xx =
Columns 1 through 6 0 1.0000e-004 2.0000e-004 3.0000e-004 4.0000e-004 5.0000e-004
Columns 7 through 11
6.0000e-004 7.0000e-004 8.0000e-004 9.0000e-004 1.0000e-003 yy =
Columns 1 through 6
1.0000e+000 9.0484e-001 8.1873e-001 7.4082e-001 6.7032e-001 6.0653e-001
Columns 7 through 11
5.4881e-001 4.9659e-001 4.4933e-001 4.0657e-001 3.6788e-001 d =
Columns 1 through 6
1.0000e+000 9.0484e-001 8.1873e-001 7.4081e-001 6.7031e-001 6.0651e-001
Columns 7 through 11
5.4878e-001 4.9654e-001 4.4927e-001 4.0649e-001 3.6778e-001
（2）预测-矫正Euler算法 程序：
function[x,y]=Euler_r(ydot_fun,x0,y0,h,N)x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;for n=1:N x(n+1)=x(n)+h;ybar=y(n)+ h*feval(ydot_fun,x(n),y(n))y(n+1)=y(n)+h/2*feval(ydot_fun,x(n),y(n))+ feval(ydot_fun,x(n+1),ybar)end
format short e ydot_fun=inline('-1000*(y-x*x)+2*x','x','y')[x,y]=Euler_r(ydot_fun,0,0,0.001,10)xx=0:0.0001:0.001 yy=exp(-1000*xx)+xx.*xx d=abs(y-yy)plot(xx,d)title('Eular矫正方法误差趋势图 h=0.001')运行结果： h=0.1 x =
Columns 1 through 6 0 1.0000e-001 2.0000e-001 3.0000e-001 4.0000e-001 5.0000e-001
Columns 7 through 11
6.0000e-001 7.0000e-001 8.0000e-001 9.0000e-001 1.0000e+000 y =
Columns 1 through 6 0 1.0200e+001 1.0083e+006 9.9774e+010 9.8728e+015 9.7692e+020
Columns 7 through 11
9.6667e+025 9.5653e+030 9.4650e+035 9.3657e+040 9.2675e+045 xx =
Columns 1 through 6 0 1.0000e-001 2.0000e-001 3.0000e-001 4.0000e-001 5.0000e-001
Columns 7 through 11
6.0000e-001 7.0000e-001 8.0000e-001 9.0000e-001 1.0000e+000 yy =
Columns 1 through 6
1.0000e+000 1.0000e-002 4.0000e-002 9.0000e-002 1.6000e-001 2.5000e-001
Columns 7 through 11
3.6000e-001 4.9000e-001 6.4000e-001 8.1000e-001 1.0000e+000 d =
Columns 1 through 6
1.0000e+000 1.0190e+001 1.0083e+006 9.9774e+010 9.8728e+015 9.7692e+020
Columns 7 through 11
9.6667e+025 9.5653e+030 9.4650e+035 9.3657e+040 9.2675e+045
h=0.01 x =
Columns 1 through 6 0 1.0000e-002 2.0000e-002 3.0000e-002 4.0000e-002 5.0000e-002
Columns 7 through 11
6.0000e-002 7.0000e-002 8.0000e-002 9.0000e-002 1.0000e-001 y =
Columns 1 through 6 0 1.2000e-001 1.0788e+003 9.7045e+006 8.7302e+010 7.8537e+014
Columns 7 through 11
7.0652e+018 6.3558e+022 5.7177e+026 5.1436e+030 4.6272e+034 xx =
Columns 1 through 6 0 1.0000e-002 2.0000e-002 3.0000e-002 4.0000e-002 5.0000e-002
Columns 7 through 11
6.0000e-002 7.0000e-002 8.0000e-002 9.0000e-002 1.0000e-001 yy =
Columns 1 through 6
1.0000e+000 1.4540e-004 4.0000e-004 9.0000e-004 1.6000e-003 2.5000e-003
Columns 7 through 11
3.6000e-003 4.9000e-003 6.4000e-003 8.1000e-003 1.0000e-002 d =
Columns 1 through 6
1.0000e+000 1.1985e-001 1.0788e+003 9.7045e+006 8.7302e+010 7.8537e+014
Columns 7 through 11
7.0652e+018 6.3558e+022 5.7177e+026 5.1436e+030 4.6272e+034
h=0.001 x =
Columns 1 through 6 0 1.0000e-003 2.0000e-003 3.0000e-003 4.0000e-003 5.0000e-003
Columns 7 through 11
6.0000e-003 7.0000e-003 8.0000e-003 9.0000e-003 1.0000e-002
y =
Columns 1 through 6 0 3.0000e-003 6.5015e-003 1.0255e-002 1.4135e-002 1.8079e-002
Columns 7 through 11
2.2057e-002 2.6053e-002 3.0058e-002 3.4069e-002 3.8084e-002 xx =
Columns 1 through 6 0 1.0000e-003 2.0000e-003 3.0000e-003 4.0000e-003 5.0000e-003
Columns 7 through 11
6.0000e-003 7.0000e-003 8.0000e-003 9.0000e-003 1.0000e-002 yy =
Columns 1 through 6
1.0000e+000 3.6788e-001 1.3534e-001 4.9796e-002 1.8332e-002 6.7629e-003
Columns 7 through 11
2.5148e-003 9.6088e-004 3.9946e-004 2.0441e-004 1.4540e-004 d =
Columns 1 through 6
1.0000e+000 3.6488e-001 1.2884e-001 3.9541e-002 4.1968e-003 1.1316e-002
Columns 7 through 11
1.9542e-002 2.5092e-002 2.9658e-002 3.3864e-002 3.7939e-002
h=0.0001 x =
Columns 1 through 6 0 1.0000e-004 2.0000e-004 3.0000e-004 4.0000e-004 5.0000e-004
Columns 7 through 11
6.0000e-004 7.0000e-004 8.0000e-004 9.0000e-004 1.0000e-003 y =
Columns 1 through 6 0 2.1000e-004-1.8838e-001 1.6937e+002-1.5227e+005 1.3690e+008
Columns 7 through 11-1.2308e+011 1.1065e+014-9.9482e+016 8.9439e+019-8.0410e+022 xx =
Columns 1 through 6 0 1.0000e-004 2.0000e-004 3.0000e-004 4.0000e-004 5.0000e-004
Columns 7 through 11
6.0000e-004 7.0000e-004 8.0000e-004 9.0000e-004 1.0000e-003 yy =
Columns 1 through 6
1.0000e+000 9.0484e-001 8.1873e-001 7.4082e-001 6.7032e-001 6.0653e-001
Columns 7 through 11
5.4881e-001 4.9659e-001 4.4933e-001 4.0657e-001 3.6788e-001 d =
Columns 1 through 6
1.0000e+000 9.0463e-001 1.0071e+000 1.6862e+002 1.5227e+005 1.3690e+008
Columns 7 through 11
1.2308e+011 1.1065e+014 9.9482e+016 8.9439e+019 8.0410e+022
现象和结论
1.对于两种方法，计算步长对计算结果稳定性有较大影响，步长选取不当时计算误差会很大；
2.在一定范围内，计算步长越小，计算结果越稳定，结果误差越小；
3.但是并不是计算步长越小越好，计算步长取值越小计算次数就会增大，初始数据误差积累越严重，也会导致结果失真。
