数值分析作业MATLABWord格式文档下载.docx
《数值分析作业MATLABWord格式文档下载.docx》由会员分享，可在线阅读，更多相关《数值分析作业MATLABWord格式文档下载.docx（12页珍藏版）》请在冰点文库上搜索。
answer1
34
3.00003.5000
3.00003.2500
3.12503.2500
3.12503.1875
3.12503.1563
3.14063.1563
3.14063.1484
3.14453.1484
3.14453.1465
3.14553.1465
3.14603.1465
3.14603.1462
3.14613.1462
3.14623.1462
最终结果为:
3.1462
2.试编写MATLAB函数实现Newton插值，要求能输出插值多项式。
对函数
在区间[-5，5]上实现10次多项式插值。
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:
ifk~=j
t=（x0（j）-x0（k））*t;
end;
s=s+y0（j）/t;
b（i）=s;
end;
t=linspace（0,0,n+1）;
n
s=linspace（0,0,n+1）;
s（n+1-i:
n+1）=b（i+1）.*poly（x0（1:
i））;
t=t+s;
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=
-（*x^10）/12928+x^9/0616+（256*x^8）/93425-x^7/076-（013693*x^6）/3421312-（3*x^5）/0+（36624*x^4）/93425-（5*x^3）/0-（*x^2）/+（7*x）/0+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（'
原函数'
十次牛顿插值多项式'
绘制结果：
误差计数并绘制误差图：
holdoff
ey=1./（1+4.*x.^2）-y
ey=
-2.6900-3.9887-4.3403-4.0857-3.4804-2.7109-1.9077-1.1575-0.5128-0.00000.37330.6130
0.73390.75580.70100.59210.45020.29430.14010.0000-0.1169-0.2051-0.2617-0.2870
-0.2832-0.2542-0.2053-0.1424-0.0719-0.00000.06740.12540.16960.19710.20620.1962
0.16790.12340.06600.0000-0.0691-0.1349-0.1902-0.2270-0.2379-0.2171-0.1649-0.0928
-0.02710-0.0271-0.0928-0.1649-0.2171-0.2379-0.2270-0.1902-0.1349-0.06910.0000
0.06600.12340.16790.19620.20620.19710.16960.12540.06740.0000-0.0719-0.1424
-0.2053-0.2542-0.2832-0.2870-0.2617-0.2051-0.11690.00000.14010.29430.45020.5921
0.70100.75580.73390.61300.37330.0000-0.5128-1.1575-1.9077-2.7109-3.4804-4.0857
-4.3403-3.9887-2.6900
plot（x,ey）
xlabel（'
ylabel（'
ey'
title（'
Runge现象误差图'
3.应用牛顿迭代法于方程f（x）-1=0,导出平方根的迭代公式，用此公式计算.
算法：
1-115/x^2'
f1=inline（'
230/x^3'
x0=10;
er=1;
k=0;
0.00001
x=x0-f（x0）/f1（x0）;
er=abs（x-x0）
x0=x;
disp（[x]）;
求解结果：
answer9
er=0.6522
10.6522
er=0.0709
10.7231
er=7.1604e-04
10.7238
er=7.1729e-08
最终结果：
4.实验数据
使用次数x
容积y
2
106.42
11
110.59
3
108.26
12
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
选用双曲线
对数据进行拟合，使用最小二乘法求出拟合函数，做出拟合曲线图。
【解】
clear,clc;
%题目条件
x=[2356791011121416171920];
y=[106.42,108.26,109.58,109.50,109.86,110.00,109.93...
110.59,110.60,110.72,110.90,110.76,111.10,111.30];
%使用最小二乘法求出1次多项式拟合系数
a=polyfit（1./x,1./y,1）;
%绘制拟合图像
xx=0.04:
0.01:
0.5;
yy=a
（1）*xx+a
（2）;
plot（1./xx,1./yy,x,y,'
*'
holdon;
xx=-0.5:
-0.04;
plot（1./xx,1./yy）;
使用最小二乘法拟合的曲线方程为
下图为绘制出的拟合曲线，并同时将一直点用“*”表示到图中。
5、利用LU分解法解方程组
首先，编辑一个LU分解函数如下
function[L,U]=Lu（A）
%求解线性方程组的三角分解法
%A为方程组的系数矩阵
%L和U为分解后的下三角和上三角矩阵
[n,m]=size（A）;
ifn~=m
error（'
TherowsandcolumnsofmatrixAmustbeequal!
'
return;
%判断矩阵能否LU分解
forii=1:
fori=1:
ii
AA（i,j）=A（i,j）;
if（det（AA）==0）
ThematrixcannotbedividedbyLU!
%开始计算，先赋初值
L=eye（n）;
U=zeros（n,n）;
%计算U的第一行，L的第一列
U（1,i）=A（1,i）;
L（i,1）=A（i,1）/U（1,1）;
%计算U的第r行，L的第r列
fori=2:
forj=i:
i-1
M（k）=L（i,k）*U（k,j）;
U（i,j）=A（i,j）-sum（M）;
forj=i+1:
M（k）=L（j,k）*U（k,i）;
L（j,i）=（A（j,i）-sum（M））/U（i,i）;
然后，编辑一个通过LU分解法解线性方程组的函数如下
function[L,U,x]=Lu_x（A,d）
%三角分解法求解线性方程组，LU法解线性方程组Ax=LUx=d
%A为方程组的系数矩阵
%d为方程组的右端项
%x为线性方程组的解
AA（i,j）=A（i,j）;
[L,U]=Lu（A）;
%直接调用自定义函数，首先将矩阵分解，A=LU
%设Ly=d由于L是下三角矩阵,所以可求y（i）
y
（1）=d
（1）;
n
forj=1:
d（i）=d（i）-L（i,j）*y（j）;
y（i）=d（i）;
%设Ux=y,由于U是上三角矩阵，所以可求x（i）
x（n）=y（n）/U（n,n）;
fori=（n-1）:
-1:
1
forj=n:
i+1
y（i）=y（i）-U（i,j）*x（j）;
x（i）=y（i）/U（i,i）;
然后，n=5时，调用自定义函数
[L,U,x]=Lu_x（A,a）
解出：
x=0.99896551.0000000000326090.99616251.0000000000199840.9996114
[L,U,x]=Lu_x（B,b）
x=0.99999761.0000000000003420.99988191.0000000000014770.9999388
升级会员 