数值分析大作业三四五六七完整版Word格式文档下载.docx
《数值分析大作业三四五六七完整版Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《数值分析大作业三四五六七完整版Word格式文档下载.docx（25页珍藏版）》请在冰豆网上搜索。
flagl=O;
m=m+l;
xO=xl;
ifflagl==labs（xO）>
=ep
flag=O;
fprintf赧人的sigma值为：
%f\n*,sigma）;
2.
求下列方程的非零根
Matlab程序为:
（1）主程序
clcformatlongx0=765;
N=100;
errorlim=10"
（-5）;
x=xO-f（xO）/subs（df（）,xO）;
n=l;
whilen<
N
x=xO-f（xO）/subs（df0,xO）;
ifabs（x-xO）>
errorlimn=n+l;
else
break;
xO=x;
disp（[‘迭代次数：
n=*,num2str（n）J）
disp（['
所求IE零根:
正根xl='
num2str（x）,'
负根x2='
num2str（-x）J）
（2）子函数非线性函数f
functiony=f（x）
y=log（（513+*x）/*x））-x/（1400*;
（3）•函数非线性函数的•阶导数df
functiony=df（）
symsxl
y=log（（513+*xl）/*xl））-xl/（1400*;
y=diff（y）;
运行结果如下：
迭代次数:
n=5
所求非零根：
正根xU负根x2=
大作业四
试编写MATLAB函数实现Newton插值，要求能输
出插值多项式.对函数/*（%）=―在区间[-5,5]1+4Q
上实现10次多项式插值.
分析：
（1）输出插值多项式。
（2）在区间[-5,5]内均匀插入99个节点，计算这些节点上函数f（x）的近似值，并在同一张图上画出原函数和插值多项式的图形。
（3）观察龙格现象，计算插值函数在各节点处的误差，并画出误差图。
Matlab程序代码如下：
%此函数实现y=l/（l+4*x*2）的n次Newton插值，n由调用函数时指定
%函数输出为插值结果的系数向量（行向量〉和插值多项式
function[ty]=func5（n）
xO=linspace（-5,5,n十1）'
y0=l./（l.+4.*x0."
2）;
b=zeros（1,n+1）;
fori=l:
n+l
s=0;
forj=l:
i
t=l；
fork=l:
ifk~=j
t=（xO（j）-xO（k））*t;
end;
s=s+yO（j）/t;
b（i）=s;
t=linspace（0,0,n+1）;
n
s=linspace（0,0,n+1）;
s（n+l-i:
n+1）=b（i+1）.*poly（x0（l:
i））;
t=t+s;
t（n+l）=t（n+l）+b（l）;
y=poly2sym（t）;
10次插值运行结果：
[bY]=func5（10）
b二
Columns1through4
Columns5through8
Columns9through11
Y二
b为插值多项式系数向量，Y为插值多项式。
插值近似值：
xl=linspace（-5,5,101）;
x=xl（2:
100）;
y=polyval（b,x）
y二
Columns1through12
绘制原函数和拟合多项式的图形代码:
plot（x,1./（1+4.*x・"
2））
holdall
plot（x,y,'
r'
）
xlabel（'
X’）ylabel（'
Y'
）title（'
Runge现象'
）gtextC原函数'
）gtext十次牛顿插值多项式'
）詁差计数并绘制误差图：
holdoff
ey=l./（1+4・*x."
2）-y
ey=
Columns13through24
Columns25through36
Columns37through48
Columns49through60
Columns61through72
Columns73through84
Columns85through96
plot（x,ey）
xlabelfX'
ylabel（'
ey'
title（'
Runge现象误差图'
输出结果为：
大作业五
Mat1ab程序为：
x=[-520,-280,,-7&
,0,,,78,,280,520]'
y=[0,-30,-36,-35,,,0,,,35,36,30,0]'
n=13;
喘求解M
fori=1:
l:
n-l
h（i）=x（i+l）-x（i）;
fori=2:
a（i）=h（i-l）/（h（i-l）+h（i））;
b（i）=l-a（i）;
c（i）=6*（（y（i+l）-y（i））/h（i）-（y（i）-y（i-l））/h（i-l））/（h（i-l）+h（i））;
a（n）=h（n-l）/（h（l）+h（n-l））;
b（n）=h（l）/（h（l）+h（n-l））;
c（n）=6/（h（l）+h（n-l））*（（y⑵-y（l））/h⑴-（y（n）-y（n-l））/h（n-l））;
A（l,1）=2;
A（l,2）=b
（2）;
A（l,n-1）=a
（2）;
A（n-1,n-2）=a（n）;
A（n-1,n-1）=2;
A（n-1,1）=b（n）;
1:
n-2
A（iti）=2;
A（i,i+1）=b（i+l）;
A（i,i-1）=a（i+l）;
end
C=c（2:
n）;
m=A\C;
M（l）=m（n-l）;
M（2:
n）=m;
xx=-520:
:
520;
51
forj=1:
ifx（j）<
=xx（i）&
xx（i）<
x（j-^l）
yy（i）=M（j+l）*（xx（i）-x（j）厂3/（6*h（j））-M（j）*（xx（i）-x（j+1））*3/（6*h（j））+（y（j+l）-M（j+1）*h（j）*2/6）*（xx（i）-x（j））/h（j）-（y（j）-M（j）*h（j）"
2/6）*（xx（i）-x（j+l））/h（j）;
end;
fori=52:
101
yy（i）=-yy（102-i）;
50
xx（i）=-xx（i）;
plot（xx,yy）;
holdon;
n/2
x（i）=-x（i）;
bd‘）；
titleC机翼外形曲线'
）;
输出结果：
运行文件，得到
2.
（1）编制求第一型3次样条插值函数的通用程序；
（2）已知汽车门曲线型值点的数据如下：
（1）Matlab编制求第一型3次样条插值函数的通用程序：
function[Sx]=Threch（X,Y,dyO,dyn）
%X为输入变量x的数值
喘Y为函数值y的数值
瓯dyO为左端•阶导数值
%dyn为右端"
阶导数值
瓯Sx为输出的函数农达式
n=length（X）~l;
d=zeros（n+1,1）;
h=zeros（1,n-1）;
fl=zeros（1,n-1）;
f2=zeros（l,n-2）;
fori=l:
n喘求函数的•阶筮瀚
h（i）=X（i+1）-X（i）;
fl（i）=（Y（i+l）-Y（i））/h（i）;
n嗚求函数的二阶签商
f2（i）=（f1（i）-fl（i-1））/（X（i+1）-X（i-1））;
d（i）=6*f2（i）;
d（l）=6*（fl（l）-dy0）/h（l）;
d（n+l）=6*（dyn-fl（n-l））/h（n~l）;
%赋初值
A=zeros（n+1,n+1）;
B=zeros（1,n-1）;
C=zeros（1,n-1）;
B（i）=h（i）/（h（i）+h（i+l））;
C（i）=l-B（i）;
A（l,2）=1；
A（n+l,n）=1;
A（i,i）=2;
A（i,i-1）=B（i-l）;
A（ifi+1）=C（i-l）;
M=A\d;
symsx;
Sx（i）=collect（Y（i）+（fl（i）-（M（i）/3+M（i+l）/6）*h（i））*（x-X（i））+M（i）/2*（x-X（i））*2+（M（i+l）-M（i））/（6*h（i））*（x-X（i））*3）;
digits（4）;
Sx（i）=vpa（Sx（i））;
dispCS（x）='
%s（%d,%d）\n,char（Sx（i））,X（i）,X（i+1））;
S=zeros（1,n）;
x=X⑴十；
S（i）=Y（i）+（f1（i）-（M（i）/3+M（i+1）/6）*h（i））*（x-X（i））+M（i）/2*（x-X（i））*2+（M（i+1）-M（i））/（6*h（i））*（x-X（i）厂3;
disp（*S（i+‘）;
dispCiX（i+S（i+'
fprintfC%d%.4f%.4f\n,i,X（i）+,S（i））;
»
X=[012345678910];
Y=[];
Threch（X,Y,,
S（x）=
*x
—
*x"
2
-*x"
3+
（0,1）
（1,2）
S（x）二
（2,3）
'
一*x
（3,4）
+*x"
3-
（4,5）
（5,6）
（6,7）
（7,8）
（&
9）
（9,10）
S（i+
•
1
X（i+
4
5
6
7
8
9
10ans
[-*x"
3-*x"
2+*x+，-*x"
2+*x+,-*x"
3+*x"
2-*x+,*x"
2+*x-,-*x"
2-*x+,*x"
2+*x-,*x"
2+*x一]
大作业六
1、炼钢厂出钢时所用的圣刚睡的钢包，在使用过程中111于钢液及炉渣对包衬耐火材料的侵蚀，使其容积不断增大，经试验，钢包的容积与相应的使用次数的数据如下：
（使用次数X,容积y）
xy
f/?
*1/x
O
o
g
选用
双曲
线
l/y=a
对使用最
乘法
数据
进行
拟合。
functiona=nihehanshu（）
x0=[2356791011121416171920];
y0=[1;
A=zeros（2,2）;
B=zeros（2,1）;
a=zeros（2,1）;
x=l・/x0;
y=l./y0;
A（l,1）=14;
A（l,2）=sum（x）;
A（2,1）=A（1,2）;
A（2,2）=sum（x.*2）;
B（l）=sum（y）;
B
（2）=sum（x.*y）;
a=A\B;
y=l./（a（l）+a
（2）*1./x0）;
subplot（1,2,2）;
plot（xO,yO-y,'
bd-'
titleC拟合曲线误差'
subplot（1,2,1）;
plot（xO,yO,'
go'
x=2:
20;
y=l./（a（l）+a
（2）*1./x）;
plot（x,y,'
r*-'
legend（'
散点,，'
拟介曲线图1/y=a
（1）+a
（2）*l/x'
title（'
最小二乘法拟介曲线'
试验所求的系数为：
nihehanshu
ans=
则拟合曲线为丄=0.008973276262446+0.000841690019705-
y%
拟合曲线图、散点图、误差图如下：
2、下面给出的是乌鲁木齐最近1个月早晨7:
00左右（新媼时间）的天气预报所得到的温度，按照数据找出任意次曲线拟合方程和它的图像。
用MATLAB编程对上述数据进行最小二乘拟合。
2008年10月--11月26日
10
温度
11
12
13
14
天数
15
16
17
18
19
20
盍
21
22
23
24
25
26
27
28
29
30
Matlab的程序如下：
x=[l:
30];
y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,&
9,11,9,7,6,5,3,1];
al=polyfit（x,5%3）%三次多项式拟合％
a2=polyfit（x,y,9）%九次多项式拟合%
a3=polyfit（x,y,15）%十五次多项式拟合%
bl=polyval（al,x）
b2=polyval（a2,x）
b3=polyval（a3,x）
rl=sum（（y-bl）.'
2）滋三次多项式i吴差平方和％
r2=sum（（y-b2）.*2）%九次次多项式误差平方和％
r3=sum（（y-b3）."
2）弔十五次多项式误差平方和％
plot（x,y/**）%用*画出x,y图像%
holdon
plot（x,bl,'
r）%用红色线画出x,bl图像%
plot（x,b2,'
g）%用绿色线画出x,b2图像%
plot（x,b3,'
b:
o‘）弘用蓝色o线画出x,b3图像%
试验结果为：
al=
Columns1through2
Columnsa2=
Columns
Columnsa3=
Columnsbl=
Columnsb2=
through4
through2
through4through6through8through10
through6
through8through10
through12
through14
through16
through18
through20
through22
through24
through26
through28
through30
throughthroughthrough
b3=
through8
through10
through
rl=r2=r3=
其中的最后图像为:
大作业七
用Romberg求积法计算积分
的近似值，要求误差不超过①5x1旷
Matlab程序为：
玄被积函数m文件：
functionFx=fx（x）
Fx=l/（l+100*x*x）;
%Romberg求积法计算积分的通用程序
functionRomberg0
a=input（'
请输入积分下限:
a='
b=input（'
请输入积分I•.限：
b='
eps=input（，请输入允许持度：
eps='
%=====计算Tn=====%
functionTn=T（n）
Tn=0;
h=（b-a）/n;
x=zeros（1,n+1）;
x（k）=a+（k-l）*h;
Tn=Tn+h*（fx（x（j））+fx（x（j+1）））/2;
%=======计算Sn======%
functionSn=S（n）
Sn=4/3*T（2*n）-l/3*T（n）;
%========计算Cn=====%
functionCn=C（n）
Cn=16/15*S（2*n）-1/15*S（n）;
%=====计算Rn=====%
functionRn=R（n）
Rn=64/63*C（2*n）-1/63*C（n）;
賀====d十算满足允许精度的Rn，并打印输｛!
«
========%
i=l；
flag=l；
ifabs（R（2"
i）-R（2"
（i-l）））/255<
eps
flag=0
升级会员 