数值分析大作业三四五六七.docx
《数值分析大作业三四五六七.docx》由会员分享，可在线阅读，更多相关《数值分析大作业三四五六七.docx（33页珍藏版）》请在冰点文库上搜索。
数值分析大作业三四五六七
大作业三
1.给定初值
及容许误差
，编制牛顿法解方程f（x）=0的通用程序.
解：
Matlab程序如下：
函数m文件：
functionFu=fu（x）
Fu=x^3/3-x;
end
函数m文件：
functionFu=dfu（x）
Fu=x^2-1;
end
用Newton法求根的通用程序
clear;
x0=input（'请输入初值x0:
'）;
ep=input（'请输入容许误差:
'）;
flag=1;
whileflag==1
x1=x0-fu（x0）/dfu（x0）;
ifabs（x1-x0）flag=0;
end
x0=x1;
end
fprintf（'方程的一个近似解为：
%f\n',x0）;
寻找最大δ值的程序：
clear
eps=input（'请输入搜索精度：
'）;
ep=input（'请输入容许误差：
'）;
flag=1;
k=0;
x0=0;
whileflag==1
sigma=k*eps;
x0=sigma;
k=k+1;
m=0;
flag1=1;
whileflag1==1&&m<=10^3
x1=x0-fu（x0）/dfu（x0）;
ifabs（x1-x0）flag1=0;
end
m=m+1;
x0=x1;
end
ifflag1==1||abs（x0）>=ep
flag=0;
end
end
fprintf（'最大的sigma值为：
%f\n',sigma）;
2．求下列方程的非零根
解：
Matlab程序为：
（1）主程序
clear
clc
formatlong
x0=765;
N=100;
errorlim=10^（-5）;
x=x0-f（x0）/subs（df（）,x0）;
n=1;
whilenx=x0-f（x0）/subs（df（）,x0）;
ifabs（x-x0）>errorlim
n=n+1;
else
break;
end
x0=x;
end
disp（['迭代次数:
n=',num2str（n）]）
disp（['所求非零根:
正根x1=',num2str（x）,'负根x2=',num2str（-x）]）
（2）子函数非线性函数f
functiony=f（x）
y=log（（513+*x）/*x））-x/（1400*;
end
（3）子函数非线性函数的一阶导数df
functiony=df（）
symsx1
y=log（（513+*x1）/*x1））-x1/（1400*;
y=diff（y）;
end
运行结果如下：
迭代次数:
n=5
所求非零根:
正根x1=负根x2=
大作业四
分析：
（1）输出插值多项式。
（2）在区间[-5,5]内均匀插入99个节点，计算这些节点上函数f（x）的近似值，并在同一张图上画出原函数和插值多项式的图形。
（3）观察龙格现象，计算插值函数在各节点处的误差，并画出误差图。
解：
Matlab程序代码如下：
%此函数实现y=1/（1+4*x^2）的n次Newton插值，n由调用函数时指定
%函数输出为插值结果的系数向量（行向量）和插值多项式
function[ty]=func5（n）
x0=linspace（-5,5,n+1）';
y0=1./（1.+4.*x0.^2）;
b=zeros（1,n+1）;
fori=1:
n+1
s=0;
forj=1:
i
t=1;
fork=1:
i
ifk~=j
t=（x0（j）-x0（k））*t;
end;
end;
s=s+y0（j）/t;
end;
b（i）=s;
end;
t=linspace（0,0,n+1）;
fori=1:
n
s=linspace（0,0,n+1）;
s（n+1-i:
n+1）=b（i+1）.*poly（x0（1:
i））;
t=t+s;
end;
t（n+1）=t（n+1）+b
（1）;
y=poly2sym（t）;
10次插值运行结果：
[bY]=func5（10）
b=
Columns1through4
Columns5through8
Columns9through11
Y=
-（735*x^10）/+x^9/9551616+（256*x^8）/93425-x^7/846976-（013693*x^6）/3421312-（3*x^5）/0+（36624*x^4）/93425-（5*x^3）/0-（511*x^2）/+（7*x）/0+1
b为插值多项式系数向量，Y为插值多项式。
插值近似值：
x1=linspace（-5,5,101）;
x=x1（2:
100）;
y=polyval（b,x）
y=
Columns1through12
Columns13through24
Columns25through36
Columns37through48
Columns49through60
Columns61through72
Columns73through84
Columns85through96
Columns97through99
绘制原函数和拟合多项式的图形代码：
plot（x,1./（1+4.*x.^2））
holdall
plot（x,y,'r'）
xlabel（'X'）
ylabel（'Y'）
title（'Runge现象'）
gtext（'原函数'）
gtext（'十次牛顿插值多项式'）
绘制结果：
误差计数并绘制误差图：
holdoff
ey=1./（1+4.*x.^2）-y
ey=
Columns1through12
Columns13through24
Columns25through36
Columns37through48
Columns49through60
0
Columns61through72
Columns73through84
Columns85through96
Columns97through99
plot（x,ey）
xlabel（'X'）
ylabel（'ey'）
title（'Runge现象误差图'）
输出结果为：
大作业五
解：
Matlab程序为：
x=[-520,-280,,-78,,,0,,,78,,280,520]';
y=[0,-30,-36,-35,,,0,,,35,36,30,0]';
n=13;
%求解M
fori=1:
1:
n-1
h（i）=x（i+1）-x（i）;
end
fori=2:
1:
n-1
a（i）=h（i-1）/（h（i-1）+h（i））;
b（i）=1-a（i）;
c（i）=6*（（y（i+1）-y（i））/h（i）-（y（i）-y（i-1））/h（i-1））/（h（i-1）+h（i））;
end
a（n）=h（n-1）/（h
（1）+h（n-1））;
b（n）=h
（1）/（h
（1）+h（n-1））;
c（n）=6/（h
（1）+h（n-1））*（（y
（2）-y
（1））/h
（1）-（y（n）-y（n-1））/h（n-1））;
A（1,1）=2;
A（1,2）=b
（2）;
A（1,n-1）=a
（2）;
A（n-1,n-2）=a（n）;
A（n-1,n-1）=2;
A（n-1,1）=b（n）;
fori=2:
1:
n-2
A（i,i）=2;
A（i,i+1）=b（i+1）;
A（i,i-1）=a（i+1）;
end
C=c（2:
n）;
C=C';
m=A\C;
M
（1）=m（n-1）;
M（2:
n）=m;
xx=-520:
:
520;
fori=1:
51
forj=1:
1:
n-1
ifx（j）<=xx（i）&&xx（i）break;
end
end
yy（i）=M（j+1）*（xx（i）-x（j））^3/（6*h（j））-M（j）*（xx（i）-x（j+1））^3/（6*h（j））+（y（j+1）-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+1））/h（j）;
end;
fori=52:
101
yy（i）=-yy（102-i）;
end;
fori=1:
50
xx（i）=-xx（i）;
end
plot（xx,yy）;
holdon;
fori=1:
1:
n/2
x（i）=-x（i）;
end
plot（x,y,'bd'）;
title（'机翼外形曲线'）;
输出结果：
运行文件，得到
2.
（1）编制求第一型3次样条插值函数的通用程序；
（2）已知汽车门曲线型值点的数据如下：
解：
（1）Matlab编制求第一型3次样条插值函数的通用程序：
function[Sx]=Threch（X,Y,dy0,dyn）
%X为输入变量x的数值
%Y为函数值y的数值
%dy0为左端一阶导数值
%dyn为右端一阶导数值
%Sx为输出的函数表达式
n=length（X）-1;
d=zeros（n+1,1）;
h=zeros（1,n-1）;
f1=zeros（1,n-1）;
f2=zeros（1,n-2）;
fori=1:
n%求函数的一阶差商
h（i）=X（i+1）-X（i）;
f1（i）=（Y（i+1）-Y（i））/h（i）;
end
fori=2:
n%求函数的二阶差商
f2（i）=（f1（i）-f1（i-1））/（X（i+1）-X（i-1））;
d（i）=6*f2（i）;
end
d
（1）=6*（f1
（1）-dy0）/h
（1）;
d（n+1）=6*（dyn-f1（n-1））/h（n-1）;
%赋初值
A=zeros（n+1,n+1）;
B=zeros（1,n-1）;
C=zeros（1,n-1）;
fori=1:
n-1
B（i）=h（i）/（h（i）+h（i+1））;
C（i）=1-B（i）;
end
A（1,2）=1;
A（n+1,n）=1;
fori=1:
n+1
A（i,i）=2;
end
fori=2:
n
A（i,i-1）=B（i-1）;
A（i,i+1）=C（i-1）;
end
M=A\d;
symsx;
fori=1:
n
Sx（i）=collect（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）;
digits（4）;
Sx（i）=vpa（Sx（i））;
end
fori=1:
n
disp（'S（x）='）;
fprintf（'%s（%d,%d）\n',char（Sx（i））,X（i）,X（i+1））;
end
S=zeros（1,n）;
fori=1:
n
x=X（i）+;
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;
end
disp（'S（i+'）;
disp（'iX（i+S（i+'）;
fori=1:
n
fprintf（'%d%.4f%.4f\n',i,X（i）+,S（i））;
end
end
输出结果：
>>X=[012345678910];
>>Y=[];
>>Threch（X,Y,,
S（x）=
*x-*x^2-*x^3+（0,1）
S（x）=
*x-*x^2-*x^3+（1,2）
S（x）=
*x-*x^2-*x^3+（2,3）
S（x）=
*x^2-*x-*x^3+（3,4）
S（x）=
*x-*x^2+*x^3-（4,5）
S（x）=
*x^2-*x-*x^3+（5,6）
S（x）=
*x-*x^2+*x^3-（6,7）
S（x）=
*x^2-*x-*x^3+（7,8）
S（x）=
*x-*x^2+*x^3-（8,9）
S（x）=
*x-*x^2+*x^3-（9,10）
S（i+
iX（i+S（i+
1
2
3
4
5
6
7
8
9
10
ans=
[-*x^3-*x^2+*x+,-*x^3-*x^2+*x+,-*x^3-*x^2+*x+,-*x^3+*x^2-*x+,*x^3-*x^2+*x-,-*x^3+*x^2-*x+,*x^3-*x^2+*x-,-*x^3+*x^2-*x+,*x^3-*x^2+*x-,*x^3-*x^2+*x-]
大作业六
1、炼钢厂出钢时所用的圣刚睡的钢包，在使用过程中由于钢液及炉渣对包衬耐火材料的侵蚀，使其容积不断增大，经试验，钢包的容积与相应的使用次数的数据如下：
（使用次数x，容积y）
x
y
x
y
2
9
3
10
5
14
6
16
7
17
9
19
10
20
选用双曲线
对使用最小二乘法数据进行拟合。
解：
Matlab程序如下：
functiona=nihehanshu（）
x0=[2356791011121416171920];
y0=[];
A=zeros（2,2）;
B=zeros（2,1）;
a=zeros（2,1）;
x=1./x0;
y=1./y0;
A（1,1）=14;
A（1,2）=sum（x）;
A（2,1）=A（1,2）;
A（2,2）=sum（x.^2）;
B
（1）=sum（y）;
B
（2）=sum（x.*y）;
a=A\B;
y=1./（a
（1）+a
（2）*1./x0）;
subplot（1,2,2）;
plot（x0,y0-y,'bd-'）;
title（'拟合曲线误差'）;
subplot（1,2,1）;
plot（x0,y0,'go'）;
holdon;
x=2:
:
20;
y=1./（a
（1）+a
（2）*1./x）;
plot（x,y,'r*-'）;
legend（'散点','拟合曲线图1/y=a
（1）+a
（2）*1/x'）;
title（'最小二乘法拟合曲线'）;
试验所求的系数为：
nihehanshu
ans=
则拟合曲线为
拟合曲线图、散点图、误差图如下：
2、下面给出的是乌鲁木齐最近1个月早晨7:
00左右（新疆时间）的天气预报所得到的温度，按照数据找出任意次曲线拟合方程和它的图像。
用MATLAB编程对上述数据进行最小二乘拟合。
2008年10月--11月26日
天数
1
2
3
4
5
6
7
8
9
10
温度
9
10
11
12
13
14
13
12
11
9
天数
11
12
13
14
15
16
17
18
19
20
温度
10
11
12
13
14
12
11
10
9
8
天数
21
22
23
24
25
26
27
28
29
30
温度
7
8
9
11
9
7
6
5
3
1
解：
Matlab的程序如下：
x=[1:
1:
30];
y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,8,9,11,9,7,6,5,3,1];
a1=polyfit（x,y,3）%三次多项式拟合%
a2=polyfit（x,y,9）%九次多项式拟合%
a3=polyfit（x,y,15）%十五次多项式拟合%
b1=polyval（a1,x）
b2=polyval（a2,x）
b3=polyval（a3,x）
r1=sum（（y-b1）.^2）%三次多项式误差平方和%
r2=sum（（y-b2）.^2）%九次次多项式误差平方和%
r3=sum（（y-b3）.^2）%十五次多项式误差平方和%
plot（x,y,'*'）%用*画出x,y图像%
holdon
plot（x,b1,'r'）%用红色线画出x,b1图像%
holdon
plot（x,b2,'g'）%用绿色线画出x,b2图像%
holdon
plot（x,b3,'b:
o'）%用蓝色o线画出x,b3图像%
试验结果为：
a1=
Columns1through2
Columns3through4
a2=
Columns1through2
Columns3through4
Columns5through6
Columns7through8
Columns9through10
a3=
Columns1through2
Columns3through4
Columns5through6
Columns7through8
Columns9through10
Columns11through12
Columns13through14
Columns15through16
b1=
Columns1through2
Columns3through4
Columns5through6
Columns7through8
Columns9through10
Columns11through12
Columns13through14
Columns15through16
Columns17through18
Columns19through20
Columns21through22
Columns23through24
Columns25through26
Columns27through28
Columns29through30
b2=
Columns1through2
Columns3through4
Columns5through6
Columns7through8
Columns9through10
Columns11through12
Columns13through14
Columns15through16
Columns17through18
Columns19through20
Columns21through22
Columns23through24
Columns25through26
Columns27through28
Columns29through30
b3=
Columns1through2
Columns3through4
Columns5through6
Columns7through8
Columns9through10
Columns11through12
Columns13through14
Columns15through16
Columns17through18
Columns19through20
Columns21through22
Columns23through24
Columns25through26
Columns27through28
Columns29through30
r1=
r2=
r3=
其中的最后图像为：
大作业七
用Romberg求积法计算积分
的近似值，要求误差不超过
。
解：
Matlab程序为：
%被积函数m文件：
functionFx=fx（x）
Fx=1/（1+100*x*x）;
end
%Romberg求积法计算积分的通用程序
functionRomberg（）
clear;
a=input（'请输入积分下限：
a='）;
b=input（'请输入积分上限：
b='）;
eps=input（'请输入允许精度：
eps='）;
%========计算Tn========%
functionTn=T（n）
Tn=0;
h=（b-a）/n;
x=zeros（1,n+1）;
fork=1:
n+1
x（k）=a+（k-1）*h;
end
forj=1:
n
Tn=Tn+h*（fx（x（j））+fx（x（j+1）））/2;
end
end
%========计算Sn========%
functionSn=S（n）
Sn=4/3*T（2*n）-1/3*T（n）;
end
%========计算Cn========%
functionCn=C（n）
Cn=16/15*S（2*n）-
升级会员 