x=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+0.6651*x）/（513-0.6651*x））-x/（1400*0.0918）;
end
（3）子函数非线性函数的一阶导数df
functiony=df（）
symsx1
y=log（（513+0.6651*x1）/（513-0.6651*x1））-x1/（1400*0.0918）;
y=diff（y）;
end
运行结果如下：
迭代次数:
n=5
所求非零根:
正根x1=767.3861负根x2=-767.3861
大作业四
分析：
（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
-0.00000.00000.0027-0.0000
Columns5through8
-0.0514-0.00000.3920-0.0000
Columns9through11
-1.14330.00001.0000
Y=
-（10035*x^10）/2928+x^9/616+（256*x^8）/93425-x^7/76-（693*x^6）/1312-（3*x^5）/+（36624*x^4）/93425-（5*x^3）/-（32311*x^2）/70496+（7*x）/+1
b为插值多项式系数向量，Y为插值多项式。
插值近似值：
x1=linspace（-5,5,101）;
x=x1（2:
100）;
y=polyval（b,x）
y=
Columns1through12
2.70033.99944.35154.09743.49262.72371.92111.17150.52740.0154-0.3571-0.5960
Columns13through24
-0.7159-0.7368-0.6810-0.5709-0.4278-0.2704-0.11470.02700.14580.23600.29490.3227
Columns25through36
0.32170.29580.25040.19150.12550.0588-0.0027-0.0537-0.0900-0.1082-0.1062-0.0830
Columns37through48
-0.03900.02450.10520.20000.30500.41580.52800.63690.73790.82690.90020.9549
Columns49through60
0.98861.00000.98860.95490.90020.82690.73790.63690.52800.41580.30500.2000
Columns61through72
0.10520.0245-0.0390-0.0830-0.1062-0.1082-0.0900-0.0537-0.00270.05880.12550.1915
Columns73through84
0.25040.29580.32170.32270.29490.23600.14580.0270-0.1147-0.2704-0.4278-0.5709
Columns85through96
-0.6810-0.7368-0.7159-0.5960-0.35710.01540.52741.17151.92112.72373.49264.0974
Columns97through99
4.35153.99942.7003
绘制原函数和拟合多项式的图形代码：
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
-2.6900-3.9887-4.3403-4.0857-3.4804-2.7109-1.9077-1.1575-0.5128-0.00000.37330.6130
Columns13through24
0.73390.75580.70100.59210.45020.29430.14010.0000-0.1169-0.2051-0.2617-0.2870
Columns25through36
-0.2832-0.2542-0.2053-0.1424-0.0719-0.00000.06740.12540.16960.19710.20620.1962
Columns37through48
0.16790.12340.06600.0000-0.0691-0.1349-0.1902-0.2270-0.2379-0.2171-0.1649-0.0928
Columns49through60
-0.02710-0.0271-0.0928-0.1649-0.2171-0.2379-0.2270-0.1902-0.1349-0.06910.0000
Columns61through72
0.06600.12340.16790.19620.20620.19710.16960.12540.06740.0000-0.0719-0.1424
Columns73through84
-0.2053-0.2542-0.2832-0.2870-0.2617-0.2051-0.11690.00000.14010.29430.45020.5921
Columns85through96
0.70100.75580.73390.61300.37330.0000-0.5128-1.1575-1.9077-2.7109-3.4804-4.0857
Columns97through99
-4.3403-3.9887-2.6900
plot（x,ey）
xlabel（'X'）
ylabel（'ey'）
title（'Runge现象误差图'）
输出结果为：
大作业五
解：
Matlab程序为：
x=[-520,-280,-156.6,-78,-39.62,-3.1,0,3.1,39.62,78,156.6,280,520]';
y=[0,-30,-36,-35,-28.44,-9.4,0,9.4,28.44,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:
10.4:
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（'机翼外形曲线'）;
输出结果：
运行jywx.m文件，得到
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）+0.5;
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+0.5）'）;
disp（'iX（i+0.5）S（i+0.5）'）;
fori=1:
n
fprintf（'%d%.4f%.4f\n',i,X（i）+0.5,S（i））;
end
end
输出结果：
>>X=[012345678910];
>>Y=[2.513.304.044.705.225.545.785.405.575.705.80];
>>Threch（X,Y,0.8,0.2）
S（x）=
0.8*x-0.001486*x^2-0.008514*x^3+2.51（0,1）
S（x）=
0.8122*x-0.01365*x^2-0.004458*x^3+2.506（1,2）
S（x）=
0.8218*x-0.01849*x^2-0.003652*x^3+2.499（2,3）
S（x）=
0.317*x^2-0.1847*x-0.04093*x^3+3.506（3,4）
S（x）=
6.934*x-1.463*x^2+0.1074*x^3-5.986（4,5）
S（x）=
4.177*x^2-21.26*x-0.2686*x^3+41.01（5,6）
S（x）=
53.86*x-8.344*x^2+0.427*x^3-109.2（6,7）
S（x）=
6.282*x^2-48.52*x-0.2694*x^3+129.6（7,8）
S（x）=
14.88*x-1.643*x^2+0.06076*x^3-39.42（8,9）
S（x）=
8.966*x-0.986*x^2+0.03641*x^3-21.67（9,10）
S（i+0.5）
iX（i+0.5）S（i+0.5）
10.50002.9086
21.50003.6784
32.50004.3815
43.50004.9882
54.50005.3833
65.50005.7237
76.50005.5943
87.50005.4302
98.50005.6585
109.50005.7371
ans=
[-0.008514*x^3-0.001486*x^2+0.8*x+2.51,-0.004458*x^3-0.01365*x^2+0.8122*x+2.506,-0.003652*x^3-0.01849*x^2+0.8218*x+2.499,-0.04093*x^3+0.317*x^2-0.1847*x+3.506,0.1074*x^3-1.463*x^2+6.934*x-5.986,-0.2686*x^3+4.177*x^2-21.26*x+41.01,0.427*x^3-8.344*x^2+53.86*x-109.2,-0.2694*x^3+6.282*x^2-48.52*x+129.6,0.06076*x^3-1.643*x^2+14.88*x-39.42,0.03641*x^3-0.986*x^2+8.966*x-21.67]
大作业六
1、炼钢厂出钢时所用的圣刚睡的钢包，在使用过程中由于钢液及炉渣对包衬耐火材料的侵蚀，使其容积不断增大，经试验，钢包的容积与相应的使用次数的数据如下：
（使用次数x，容积y）
x
y
x
y
2
106.42
9
110.59
3
108.26
10
110.60
5
109.58
14
110.72
6
109.50
16
110.90
7
109.86
17
110.76
9
110.00
19
111.10
10
109.93
20
111.30
选用双曲线
对使用最小二乘法数据进行拟合。
解：
Matlab程序如下：
functiona=nihehanshu（）
x0=[2356791011121416171920];
y0=[106.42108.26109.58109.50109.86110.00109.93110.59110.60110.72110.90110.76111.10111.30];
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:
0.5:
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=
0.2446
0.9705
则拟合曲线为
拟合曲线图、散点图、误差图如下：
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）%九次多项式拟合
升级会员 