flag='failure'
return;
end
fork=i+1:
n
m=0;
forq=1:
i-1
m=m+L（k,q）*U（q,i）;
end
L（k,i）=（A（k,i）-m）/U（i,i）;
end
end
L
U
end
4.结果
>>LU_decom（A）
m=
4
n=
4
L=
1000
2100
1210
3321
U=
2426
0123
0036
0001
5.拓展
在编写程序过程中由于角标较多因此在运行过程中出现了不少角标不对的错误
题目；给出函数f（x）=1/（1+25x^2），求f（x在[-1,1]上取5个和9个等距节点，做最小二乘拟合，得出均方误差方误差。
五个节点时，matlab编码为：
首先建立M文件，并保存
functiony=f（x）
y=1/（1+25*x^2）;
end
x=[-1-0.500.51];
fori=1:
5
y（i）=f（x（i））;
end
a=polyfit（x,y,3）
symsx
f1=a
（1）*x^3+a
（2）*x^2+a（3）*x+a（4）
x=[-1-0.500.51];
fori=1:
5
E（i）=（f（x（i））-（a
（1）*x（i）^3+a
（2）*x（i）^2+a（3）*x（i）+a（4）））^2;
end
sum（E）
输出结果为
a=
-0.0000-0.6063-0.00000.5737
f1=
-4878849915647781/1298074214633706907132624082305024*x^3-1600/2639*x^2-5348847520430703/649037107316853453566312041152512*x+1514/2639
（拟合的多项式）
ans=
0.3534（均方误差）
九个点的时候，matlab编码为：
x=[-1-0.75-0.5-0.2500.250.50.751];
fori=1:
9
y（i）=f（x（i））;
end
a=polyfit（x,y,3）
symsx
f2=a
（1）*x^3+a
（2）*x^2+a（3）*x+a（4）
x=[-1-0.75-0.5-0.2500.250.50.751];
fori=1:
5
E1（i）=（f（x（i））-（a
（1）*x（i）^3+a
（2）*x（i）^2+a（3）*x（i）+a（4）））^2;
end
sum（E1）
输出结果为：
a=
-0.0000-0.56090.00000.4855
f2=
-728732707776039/2535301200456458802993406410752*x^3-1263030908712921/2251799813685248*x^2+4915246442354361/20282409603651670423947251286016*x+1093229300764671/2251799813685248
（最小二乘拟合多项式）
ans=
0.3350
（均方误差）
用复合梯形公式求积分
的值。
functioni=combinetraprl（f,a,b,eps）
%复化梯形公式求函数f在区间[a，b]上的定积分
%函数名：
f
%积分下限：
a
%积分上限：
b
%积分精度：
eps
%积分值：
i
%积分划分的子区间个数：
step
n=1;
h=（b-a）/2;
i1=0;
i2=（subs（sym（f）,findsym（sym（f））,a）+subs（sym（f）,findsym（sym（f））,b））/h;
whileabs（i2-i1）>eps
n=n+1;
h=（b-a）/n;
i1=i2;
i2=0;
fori=0:
n-1
x=a+h*i;
x1=x+h;
i2=i2+（h/2）*（subs（sym（f）,findsym（sym（f））,x）+...
subs（sym（f）,findsym（sym（f））,x1））;
end
end
四阶龙格-库塔法分别求解下列初值问题；
function[x,y]=runge_kutta（fun,x0,xt,y0,pointnum,varargin）
ifnargin<3
y0=0;
end
y（1,:
）=y0（:
）';
h=（xt-x0）/（pointnum-1）;
x=x0+[0:
pointnum]'*h;
fork=1:
pointnum
f1=h*feval（fun,x（k）,y（k,:
）,varargin{:
}）;
f1=f1（:
）';
f2=h*feval（fun,x（k）+h/2,y（k,:
）+f1/2,varargin{:
}）;
f2=f2（:
）';
f3=h*feval（fun,x（k）+h/2,y（k,:
）+f2/2,varargin{:
}）;
f3=f3（:
）';
f4=h*feval（fun,x（k）+h,y（k,:
）+f3,varargin{:
}）;
f4=f4（:
）';
y（k+1,:
）=y（k,:
）+（f1+2*（f2+f3）+f4）/6;
end
function[x,y]=euler（fun,x0,xt,y0,pointnum）
ifnargin<4
y0=0;
end
h=（xt-x0）/pointnum;
x=x0+[0:
pointnum]'*h;
y（1,:
）=y0（:
）';
fork=1:
pointnum
f=feval（fun,x（k）,y（k,:
））;
f=f（:
）';
y（k+1,:
）=y（k,:
）+h*f;
end
Runge_kutta方法
>>x0=0;
xt=1;
fun=inline（'9*y/（1+2*x）','x','y'）;
y0=1;
pointnum=10;
[x1,yr]=runge_kutta（fun,x0,xt,y0,pointnum）
一．运行结果：
1）Runge_kutta方法：
x1=
0
0.1111
0.2222
0.3333
0.4444
0.5556
0.6667
0.7778
0.8889
1.0000
1.1111
yr=
1.0000
2.4611
5.2136
9.9212
17.4198
28.7302
45.0700
67.8650
98.7600
139.6286
192.5834
小行星轨道问题：
一天文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立以太阳为原点的直角坐标系，在五个不同的对小行星作了五次观察，测得轨道上五个点的坐标数据（单位：
万公里）如下表所示：
P1
P2
P3
P4
P5
X坐标
53605
58460
62859
66662
68894
Y坐标
6026
11179
16954
23492
68894
由开普勒第一定律知，小行星轨道为一椭圆，椭圆的一般方程可表示为：
a1x2+2a2xy+a3y2+2a4x+2a5y+1=0
现需要建立椭圆的方程以供研究。
（1）分别将五个点的数据代入椭圆一般方程中，写出五个待定系数满足的等式，整理后写出线性方程组
AX=R
以及方程组的系数矩阵和右端项b；
（2）用MARLAB求低阶方程的指令A＼b求出待定系数a1,a2,a3,a4,a5.
答：
1）
2）
一．源程序：
>>A=[53605^22*53605*60266026^22*536052*6026;
58460^22*58460*1117911179^22*584602*11179;
62859^22*62859*1695416954^22*628592*16954;
66662^22*66662*2349223492^22*666622*23492;
68894^22*68894*6889468894^22*688942*68894];
>>b=[-1,-1,-1,-1,-1];
>>x=A\b
二．运行结果：
数值分析课程设计报告
摘要（中）：
本文建立在数值分析的理论基础上，能够在Matlab环境中运行，给出了理论分析、程序清单以及计算结果。
更重要的是，还有详细的对算法的框图说明。
首先运用Romberg积分方法对给出定积分进行积分，然后对得到的结果用插值方法，分别求出Lagrange插值多项式和Newton插值多项式，再运用最小二乘法的思想求出拟合多项式，最后对这些不同类型多项式进行比较，找出它们各自的优劣。
Summary（Enlish）:
Thisessayisbasedonnumericalanalysis,andcouldbeoperatedinMatlabenvironment,includingtheroyanalysis,programandresults.What’smore,therearedetaileddiagramwhichshowshowthealgorithmworks.ItfirstusestheintegratingmethodofRomberg,whichisanimprovedtrapezoidalintegration,tosolvethegivendefiniteintegral,thenwecreateLagrange’sinterpolationpolynomialandNewton’sinterpolationpolynomial.Andaccordingtoleastsquaremethod,curvefittingpolynomialiscreated.Atthelastpartoftheessay,wecomparethesedifferentpatternsofpolynomial,foundingtheirdistinctiveadvantagesanddisadvantages.
主题词：
Romberg积分，插值方法，Langrange插值多项式，Newton插值多项式，拟合多项式。
一.问题提出：
已知椭圆的周长可以表示成s=a
（0<
<1）,取a=1,
⑴针对
从0.1到0.9（步长h=0.1）分别求出周长s;（用Romberg积分方法）
⑵对于以上数据，求出
的插值多项式；
⑶对于⑴中数据，试用最小二乘法的思想求作拟合多项式（要求是偶次），并对这些多
项式的优劣进行比较。
二．问题解决：
依据问题出现先后顺序，对三个小问进行讨论。
⑴Romberg算法是在复化梯形求积公式的基础上，应用理查逊外推构造的一种数值积分方法。
由复合梯形公式的展开定理，得到如下关系式：
T1（h）-I=aah2+a2h4+a3h6+…+amh2m+…
其中，I=
T1（h）=Tn.
利用Richardson外推定理对T1（h）进行加速，注意这里取m=1,q=
有
利用T0（
）和T0（
）可以得到
实际上T1（h）就是复化抛物线求积公式，一般的计算公式
（m=1,2,…,k=0,1,2,…）
即
.
由于Romberg求积过程是每次把区间缩小一半，所以Romberg积分方法也叫做逐次分半加速收敛法。
Robermg求积算法的计算过程如下：
⑴取k=0,h=b-a,求
=
令1→k（k记区间[a,b]的二分次数）.
⑵求梯形值T0
，即按递推公式
计算
.
⑶求加速值，按公式
逐个求出如表的第k行其余
各元素Tj（k-j）（j=1,2,…,k）.
⑷若
<
（预先给定的精度），则终止计算，并取Tk（0）≈I,否则令k+1→k转⑵继续计算。
Romberg算法的计算过程列出如下表：
k
h（步长）
T0（k）
T1（k）
T2（k）
T3（k）
0
h
T0（0）①
1
h/2
T0
（1）②
T1（0）③
2
h/22
T0
（2）④
T1
（1）⑤
T2（0）⑥
3
h/23
T0（3）⑦
T1
（2）⑧
T2
（1）⑨
T3（0）⑩
…
…
…
…
…
…
表中①~⑩表示计算顺序，k表示二分次数。
程序框图：
编写Matlab程序Romberg.m:
functionRomberg（p,k）
M=1;
a=0;
b=2*pi;
h=b-a;
err=1;
i=0;
R=zeros（4,4）;
R（1,1）=h*（feval（'f',p,a）+feval（'f',p,b））/2;
while（err>0.0001&ii=i+1;
h=h/2;
升级会员 