break;
else
x（:
1）=x（:
2）;
forj=[1:
1:
n],
ifj==1,
s1=0;
fore1=[2:
1:
m],
s1=s1+A（1,e1）*x（e1,1）;%求和
end
x（j,2）=（b（1,1）-s1）/A（1,1）;
elseifj==m,
s2=0;
fore2=[1:
1:
m-1],
s2=s2+A（m,e2）*x（e2,2）;%求和
end
x（j,2）=（b（m,1）-s2）/A（m,m）;
else
s3=0;
fore3=[1:
1:
j-1],
s3=s3+A（j,e3）*x（e3,2）;%求和
end
fore4=[j+1:
1:
m],
s3=s3+A（j,e4）*x（e4,1）;%求和
end
x（j,2）=（b（j,1）-s3）/A（j,j）;
end
end
end
end
x1=x（:
2）;
【结果分析与讨论】
在matlab命令窗口中输入：
>>A=[5-1-1-1;-110-1-1;-1-15-1;-1-1-110;];
>>b=[-4;12;8;34;];
>>x0=[1;1;1;1];
>>GS（A,b,x0,5,0.000001）
ans=
0.996552514803221
1.998599188612620
2.998420412137192
3.999357211555304
评价：
相对于简单迭代法，赛德尔迭代法收敛速度很快，此程序可以对任意维的方程组求指定的精度的迭代结果，实用性比较强。
实验二最小二乘法拟合数据
【实验目的】
熟悉最小二乘法拟合数据的方法
【实验内容】
1.了解matlab语言
2.给出数据：
-1.00
0.50
-0.50
-0.25
0.00
10.2209
0.3295
0.8826
1.4392
2.0003
0.25
0.50
0.75
1.00
2.5645
3.1334
3.7061
4.2836
用一次、二次多项式利用最小二乘法拟合这些数据，试写出正规方程组，并求出最小平方逼近多项式。
【实验所使用的仪器设备与软件平台】
计算机、matlab7.0
【实验方法与步骤】
最小二乘法拟合算法：
S1：
输入数据列向量
与
，拟合多项式的次数
；
S2：
计算
的列数为
；
S3：
置
（其中
表示矩阵
的第1行第1列，以下的
等类似）；
S4：
计算：
S5：
利用
分解法求解
，得到
次最小平方逼近多项式的系数与多项式；
S6：
在同一坐标系中绘出原数据点与次
多项式逼近曲线。
【实验结果】
实验程序：
functiony=fmintwo（a,m,x）%调用拟合多项式的函数
y=a（1,1）;
fore1=[1:
1:
m],
y=y+（a（e1+1,1））*（x.^e1）;
end
functiona=mintwo（x,y,m）%最小二乘法拟合程序
[n,p]=size（x）;
a（m+1,1）=0;
A（m+1,m+1）=0;
b（m+1,1）=0;
c（2*m+1,1）=0;
c（1,1）=n;
fork=[2:
1:
2*m+1],
s1=0;
fore1=[1:
1:
n],
s1=s1+（x（e1,1））^（k-1）;
end
c（k,1）=s1;
end
fori=[1:
1:
m+1],
forj=[1:
1:
m+1],
A（i,j）=c（i+j-1,1）;%求正规方程组的系数矩阵
end
end
A
forh=[1:
1:
m+1],
s2=0;
fore2=[1:
1:
n],
s2=s2+y（e2,1）*（（x（e2,1））^（h-1））;
end
b（h,1）=s2;%求正规方程组的
向量
end
b
[LDa]=ja2911（A,b）;%调用
分解法求解正规矩阵的程序
z=[min（x）:
0.1:
max（x）];
plot（x,y,'*',z,fmintwo（a,m,z）,'r'）;%绘出原数据表示的点和拟合曲线
附
分解法求系数矩阵为对称矩阵的方程组程序：
function[LDx]=ja2911（A,b）
[n,n]=size（A）;
L=eye（size（A））;
D=zeros（size（A））;
fori=[2:
1:
n],
A（i,1）=A（i,1）/A（1,1）;
forj=[2:
1:
n],
ifi>j;
s1=0;
fore1=[1:
1:
j-1],
s1=s1+A（e1,e1）*A（i,e1）*A（j,e1）;
end
A（i,j）=（A（i,j）-s1）/A（j,j）;
elseifi==j,
s2=0;
fore2=[1:
1:
i-1],
s2=s2+A（e2,e2）*A（i,e2）*A（i,e2）;
end
A（i,j）=A（i,j）-s2;
end
end
end
forh=[1:
1:
n],
D（h,h）=A（h,h）;
end
fork=[2:
1:
n],
fort=[1:
1:
k-1],
L（k,t）=A（k,t）;
end
end
y（n,1）=0;
x（n,1）=0;
y（1,1）=b（1,1）/D（1,1）;
fora=[2:
1:
n],
s3=0;
fore3=[1:
1:
a-1],
s3=s3+L（a,e3）*D（e3,e3）*y（e3,1）;
end
y（a,1）=（b（a,1）-s3）/D（a,a）;
end
x（n,1）=y（n,1）;
forbb=[n-1:
-1:
1],
s4=0;
fore4=[bb+1:
1:
n],
s4=s4+L（e4,bb）*x（e4,1）;
end
x（bb,1）=y（bb,1）-s4;
end
【结果分析与讨论】
在matlab命令窗口中输入：
>>xx=[-1;-0.75;-0.5;-0.25;0;0.25;0.5;0.75;1;];
>>yy=[-0.2209;0.3295;0.8826;1.4392;2.0003;2.5645;3.1334;…
3.7061;4.2836;];
>>mintwo（xx,yy,1）
A=
9.0000000000000000
03.750000000000000
b=
18.118299999999998
8.443674999999999
ans=
2.013144444444444
2.251646666666666
输出的图像：
>>mintwo（xx,yy,2）
A=
9.00000000000000003.750000000000000
03.7500000000000000
3.75000000000000002.765625000000000
b=
18.118299999999998
8.443674999999999
7.586956250000000
ans=
2.000100432900433
2.251646666666666
.0313********
输出的图像：
一次多项式拟合的正规方程组为：
一次拟合多项式为：
二次多项式拟合的正规方程组为：
二次拟合多项式为：
评价：
此程序实现了对数据利用最小二乘法进行
次多项式逼近，并且可以输出原数据表示的点和拟合曲线在同一坐标系中的图像，可以很直观的进行观察拟合的程度。
实验三龙贝格求积分
【实验目的】
熟悉龙贝格求积分的方法
【实验内容】
1.了解matlab语言
2.用龙贝格方法计算积分：
要求误差不超过
。
【实验所使用的仪器设备与软件平台】
计算机、matlab7.0
【实验方法与步骤】
龙贝格求积分算法：
S1：
输入区间左右端点
与
和精度
；
S2：
置
；
S3：
；
S4：
计算：
S5：
置
；
S6：
计算：
S7：
若
，则输出
，并且结束；
若
时，当
则转到s5；
当
则无输出，并且结束。
【实验结果】
实验程序：
functiony=fromberg（x）%可被调用的积分函数
y=2*exp（-（x^2））/（（pi）^（1/2））;
functionI=romberg（a,b,E）%龙贝格求积分程序
T（50,1）=0;
C（50,1）=0;
R（50,1）=0;
h（50,1）=0;
h（1,1）=（b-1）/2;
T（1,1）=h（1,1）*（fromberg（a）+fromberg（b））;
fork1=[2:
1:
5],
h（k1,1）=（b-a）/（2^（k1-1））;
g1=0;
fore1=[1:
1:
2^（k1-2）],
g1=g1+fromberg（a+（2*e1-1）*h（k1,1））;
end
T（k1,1）=（T（k1-1,1）/2）+h（k1,1）*g1;%计算
end
fori1=[1:
1:
4],
S（i1,1）=（4*T（i1+1,1）-T（i1,1））/3;%计算
end
fori2=[1:
1:
3],
C（i2,1）=（16*S（i2+1,1）-S（i2,1））/15;%计算
end
fori3=[1:
1:
2],
R（i3,1）=（64*C（i3+1,1）-C（i3,1））/63;%计算
end
fork=[1:
1:
40],
ifabs（R（k+1,1）-R（k,1））<=E,%误差判断
I=R（k+1,1）;
break;
else%如果不符合精度要求，则自动将区间分半继续计算
h（k+5,1）=（b-a）/（2^（k+4））;
g2=0;
fore2=[1:
1:
2^（k+3）],
g2=g2+fromberg（a+（2*e2-1）*h（k+5,1））;
end
T（k+5,1）=（T（k+4,1）/2）+h（k+5,1）*g2;
S（k+4,1）=（4*T（k+5,1）-T（k+4,1））/3;
C（k+3,1）=（16*S（k+4,1）-S（k+3,1））/15;
R（k+2,1）=（64*C（k+3,1）-C（k+2,1））/63;
end
end
【结果分析与讨论】
在matlab命令窗口中输入：
>>romberg（0,1,10^（-5））
ans=
0.842693582047394
评价：
龙贝格求积分方法精确度比较高，该程序很好的实现了龙贝格方法，但是程序的某些地方还可以简化，使程序精简。
实验四标准四阶R-K方法
【实验目的】
熟悉标准四阶R-K的方法
【实验内容】
1.了解matlab语言
2.用标准四阶R-K求解初值问题
的解在-0.9，-0.8，-0.7的近似值（按四位小数计算）。
【实验所使用的仪器设备与软件平台】
计算机、matlab7.0
【实验方法与步骤】
S1：
输入区间左右端点
与
、
的初值
、初始步长
、最大迭代次数
、允许误差
；
S2：
置
；
S3：
计算：
S4：
置
；
S4：
计算：
S5：
若
，则将
，输出
，并且结束；
若
，当
，将
，并转向s4；
当
，无输出，并且结束。
【实验结果】
实验程序：
functionfk=fRK（x,y）%微分方程的
函数
fk=（x^2）-（y^2）;
functionx=RK（a,b,ya,h）%用R-K方法求微分方程数值解的程序
[m,n]=size（[a:
h:
b]）;
x（n,2）=0;
x（:
1）=[a:
h:
b]';
x（1,2）=ya;
k（4,1）=0;
fori=[2:
1:
n],
k（1,1）=h*fRK（x（i-1,1）,x（i-1,2））;
k（2,1）=h*fRK（x（i-1,1）+h/2,x（i-1,2）+k（1,1）/2）;
k（3,1）=h*fRK（x（i-1,1）+h/2,x（i-1,2）+k（2,1）/2）;
k（4,1）=h*fRK（x（i-1,1）+h,x（i-1,2）+k（3,1））;
x（i,2）=x（i-1,2）+（k（1,1）+2*k（2,1）+2*k（3,1）+k（4,1））/6;
end
function[xh1]=WRK（a,b,ya,h,N,E）%R-K的步长自动选择程序
m（N,1）=0;
x1（N,1）=0;
fori=[1:
1:
N],
[lm（i,1）]=size（[a:
h/（2^（i-1））:
b]）;
end;
y=RK（a,b,ya,h）;
x1（1,1）=y（m（1,1）,2）;
y=RK（a,b,ya,h/2）;
x1（2,1）=y（m（2,1）,2）;
forj=[3:
1:
N],
ifabs（x1（j-1,1）-x1（j-2,1））<=E,
x=RK（a,b,ya,h/（2^（j-2）））;
h1=h/（2^（j-2））;
break;
else
y=RK（a,b,ya,h/（2^（j-1）））;
x（j,1）=y（m（3,1）,2）;
end
end
（该程序输出的是满足精度要求的区间分点的数值解和及相应步长）
【结果分析与讨论】
在matlab命令窗口中输入：
>>[xh]=WRK（-1,-0.5,0,0.1,20,0.00001）
x=
-1.0000000000000000
-0.9500000000000000.047503044462813
-0.9000000000000000.090047822422416
-0.8500000000000000.127736329909335
-0.8000000000000000.160727765024206
-0.7500000000000000.189********6077
-0.7000000000000000.213483856644815
-0.6500000000000000.233766*********
-0.6000000000000000.250369794944545
-0.5500000000000000.263601751827377
-0.5000000000000000.273776792527207
h=
0.050000000000000
所以微分方程的解在-0.9，-0.8，-0.7的近似值分别为0.09004、0.16072、0.21348
评价：
标准四阶R-K方法的精确度比较高，而且算法容易实现。
次程序实现了标准四阶R-K算法，通过多次调用函数的形式使程序得到了精简。
升级会员 