西南交大数值分析第二次大作业可以运行Word文档下载推荐.docx
- 文档编号:18214194
- 上传时间:2022-12-14
- 格式:DOCX
- 页数:25
- 大小:384.11KB
西南交大数值分析第二次大作业可以运行Word文档下载推荐.docx
《西南交大数值分析第二次大作业可以运行Word文档下载推荐.docx》由会员分享,可在线阅读,更多相关《西南交大数值分析第二次大作业可以运行Word文档下载推荐.docx(25页珍藏版)》请在冰豆网上搜索。
插值多项式'
yt=Lang(x,y)
ezplot(yt,[-1,5]);
运行结果:
插值结果:
Yt=12.0000
插值多项式:
yt=4.0*t^4-2.0*t^3+t^2-3.0*t+2.0
(2)构造arctanx在[1,6]基于等距节点的10次插值多项式
functionf=New(x,y,x0)
if(length(x)==length(y))
n=length(x);
c(1:
n)=0.0;
else
disp('
xº
Í
yÎ
¬
Ê
ý
²
»
µ
È
£
¡
'
);
return;
f=y
(1);
y1=0;
xx=linspace(x
(1),x(n),(x
(2)-x
(1)));
n-1)
n-i)
y1(j)=y(j+1)-y(j);
c(i)=y1
(1);
l=t;
for(k=1:
l=l*(t-k);
f=f+c(i)*l/factorial(i);
y=y1;
if(i==n-1)
(x0-x
(1))/(x
(2)-x
(1)));
>
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)];
yt=New(x,y)
ezplot(yt,[1,6]);
holdon
ezplot('
atan(t)'
[1,6])
gridon
插值多项式
yt=
1.34684*10^(-10)*t^10-6.61748*10^(-9)*t^9+1.28344*10^(-7)*t^8-0.00000104758*t^7-0.00000243837*t^6+0.000149168*t^5-0.00176296*t^4+0.0125826*t^3-0.0640379*t^2+0.250468*t+0.785398
2
(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);
holdon;
symsx
fx=1/(1+x^2);
ezplot(fx);
由图可知,三次样条插值多项式图像与原函数图像基本一致。
(2)取第一类边界条件,用三弯矩法编写MATLAB程序
程序
functionyi=cubic_spline(x,y,ydot,xi)
ny=length(y);
h=zeros(1,n);
lambda=ones(1,n);
mu=ones(1,n);
M=zeros(n,1);
d=zeros(n,1);
fork=2:
n
h(k)=x(k)-x(k-1);
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));
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));
fori=1:
A(i,i+1)=lambda(i);
A(i+1,i)=mu(i+1);
M=A\d;
ifx(k-1)<
=xi&
xi<
=x(k)
yi=M(k-1)/6/h(k)*(x(k)-xi)^3...
+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);
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:
yy=1./(1+xx.^2);
m=length(xx);
z=zeros(1,m);
m
z(i)=cubic_spline(x,y,ydot,xx(i));
o'
xx,yy,'
k:
xx,z,'
k-'
3
本题直接利用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)
Goodness
of
fit:
SSE:
0.1733
R-square:
0.9993
Adjusted
RMSE:
0.1112
(2)最小二乘法插值:
interpolant:
piecewise
polynomial
computed
from
p
Coefficients:
p
coefficient
structure
NaN
R-square(确定系数)
SSR:
Sum
squares
the
regression,即预测数据与原始数据均值之差的平方和,公式为:
SST:
Total
sum
squares,即原始数据和均值之差的平方和,公式为:
“确定系数”是定义为SSR和SST的比值,故
其实“确定系数”是通过数据的变化来表征一个拟合的好坏。
由上面的表达式可以知道“确定系数”的正常取值范围为[01],越接近1,表明方程的变量对y的解释能力越强,这个模型对数据拟合的也较好
4.
(1)16点复合梯形公式:
functionI=T_quad(x,y)
m=length(y);
h=(x(n)-x
(1))/(n-1);
a=[12*ones(1,n-2)1];
I=h/2*sum(a.*y);
End
(1)x=0:
0.0625:
y=exp(x)
I=T_quad(x,y)
x=
Columns1through7
00.06250.12500.18750.25000.31250.3750
Columns8through14
0.43750.50000.56250.62500.68750.75000.8125
Columns15through17
0.87500.93751.0000
y=
1.00001.06451.13311.20621.28401.36681.4550
1.54881.64871.75511.86821.98872.11702.2535
2.39892.55362.7183
I=
1.7188
(2)x=1:
y=sin(x)/x
1.00001.06251.12501.18751.25001.31251.3750
1.43751.50001.56251.62501.68751.75001.8125
1.87501.93752.0000
0.6116
0.6116
(2)8点的复合simpson公式:
function
I=S_quad(x,y)
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
I=
h/6*sum(a.*y)
1.7183
(2)x=1:
(3)三点的Gauss-Legendre积分公式
functionI=G_quad(fun,a,b,N)
h=(b-a)/N;
I=0;
fork=1:
t=[-sqrt(3/5)0sqrt(3/5)];
A=[5/98/95/9];
F=feval(fun,h/2*t+a+(k-1/2)*h);
I=I+sum(A.*F)
I=h/2*I;
(1)fun=inline('
exp(x)'
I=G_quad(fun,0,1,4)
2.2722
5.1898
8.9360
13.7463
(2)fun=inline('
sin(x)/x'
I=G_quad(fun,1,2,4)
1.5954
3.0143
4.2363
5.2479
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;
forn=1:
x(n+1)=x(n)+h;
y(n+1)=y(n)+h*feval(ydot_fun,x(n),y(n));
formatshorte
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:
yy=exp(-1000*xx)+xx.*xx
d=abs(y-yy)
plot(xx,d)
title('
Eular方法误差趋势图h=0.1'
h=0.1
Columns1through6
01.0000e-0012.0000e-0013.0000e-0014.0000e-0015.0000e-001
Columns7through11
6.0000e-0017.0000e-0018.0000e-0019.0000e-0011.0000e+000
Columns1through6
001.0200e+000-9.6940e+0019.6061e+003-9.5099e+005
9.4148e+007-9.3207e+0099.2274e+011-9.1352e+0139.0438e+015
yy=
1.0000e+0001.0000e-0024.0000e-0029.0000e-0021.6000e-0012.5000e-001
3.6000e-0014.9000e-0016.4000e-0018.1000e-0011.0000e+000
d=
1.0000e+0001.0000e-0029.8000e-0019.7030e+0019.6060e+0039.5099e+005
9.4148e+0079.3207e+0099.2274e+0119.1352e+0139.0438e+015
h=0.01
01.0000e-0022.0000e-0023.0000e-0024.0000e-0025.0000e-002
6.0000e-0027.0000e-0028.0000e-0029.0000e-0021.0000e-001
001.2000e-003-6.4000e-0036.7200e-002-5.8800e-001
5.3180e+000-4.7825e+0014.3047e+002-3.8742e+0033.4868e+004
xx=
01.0000e-0022.0000e-0023.0000e-0024.0000e-0025.0000e-002
1.0000e+0001.4540e-0044.0000e-0049.0000e-0041.6000e-0032.5000e-003
3.6000e-0034.9000e-0036.4000e-0038.1000e-0031.0000e-002
1.0000e+0001.4540e-0048.0000e-0047.3000e-0036.5600e-0025.9050e-001
5.3144e+0004.7830e+0014.3047e+0023.8742e+0033.4868e+004
4.9580e+0004.8315e+0014.2983e+0023.8750e+0033.4867e+004
h=0.001
01.0000e-0032.0000e-0033.0000e-0034.0000e-0035.0000e-003
6.0000e-0037.0000e-0038.0000e-0039.0000e-0031.0000e-002
003.0000e-0068.0000e-0061.5000e-0052.4000e-005
3.5000e-0054.8000e-0056.3000e-0058.0000e-0059.9000e-005
01.0000e-0032.0000e-0033.0000e-0034.0000e-0035.0000e-003
1.0000e+0003.6788e-0011.3534e-0014.9796e-0021.8332e-0026.7629e-003
2.5148e-0039.6088e-0043.9946e-0042.0441e-0041.4540e-004
1.0000e+0003.6788e-0011.3534e-0014.9788e-0021.8317e-0026.7389e-003
2.4798e-0039.1288e-0043.3646e-0041.2441e-0044.6400e-005
h=0.0001
01.0000e-0042.0000e-0043.0000e-0044.0000e-0045.0000e-004
6.0000e-0047.0000e-0048.0000e-0049.0000e-0041.0000e-003
1.0000e+0009.0484e-0018.1873e-0017.4082e-0016.7032e-0016.0653e-001
5.4881e-0014.9659e-0014.4933e-0014.0657e-0013.6788e-001
1.0000e+0009.0484e-0018.1873e-0017.4081e-0016.7031e-0016.0651e-001
5.4878e-0014.9654e-0014.4927e-0014.0649e-0013.6778e-001
(2)预测-矫正Euler算法
function[x,y]=Euler_r(ydot_fun,x0,y0,h,N)
y=zeros(1,N+1);
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)
[x,y]=Euler_r(ydot_fun,0,0,0.001,10)
0.0001:
0.001
Eular矫正方法误差趋势图h=0.001'
6.0000e-0017.0000e-0018.0000e-0019.0000e-
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 西南 交大 数值 分析 第二次 作业 可以 运行