数值积分matlab实验程序.docx
《数值积分matlab实验程序.docx》由会员分享，可在线阅读，更多相关《数值积分matlab实验程序.docx（25页珍藏版）》请在冰点文库上搜索。
数值积分matlab实验程序
这是三次样条插值，及龙格现象的图。
这里版权申明：
程序版权归崛起强（这是XXID）所有。
其中只有部分程序摘自《精通matlab科学计算》，但均作出了注解及修改。
请勿转载！
！
！
请勿用于商业用途！
！
！
用所编写的程序计算积分
，
，
复合梯形公式、复合辛普森公式、龙贝格公式、高斯列让德公式、高斯拉盖尔求积公式
function[I,n]=RecursionTrapezoid（f,a,b,error）
%n为节点数
ifnargin<3
disp（'Youmustinputthreevariablesatleast'）;
%quit;用quit不好，直接退出matlab，改为return,才能看到disp。
return;
elseifnargin==3
error=1.0e-4;%error是个精度eps,精度太高难以等待1.0e-10相当久
end
h=b-a;%步长
%y0=0.5*h*（fun（a）+fun（b））;%一步
%n=1;%%增加的节点数
%flag=1;
%while（flag）
%sum=0;%
%xs=a+0.5*h;
%fori=1:
n;%默认的做了步插值点。
翻倍插值
%sum=sum+fun（xs）;
%xs=xs+h;
%end
%y1=0.5*y0+0.5*h*sum;%y0是对前面结果的复用，sum是新节点之和，新步长0.5h
%ifabs（y1-y0）<=error
%y=y1;
%flag=0;
%end
%h=h/2;
%y0=y1;%转化到while时这句要相应修改
%n=2*n;%假如error退出，引起不一致性，n并未翻倍
%end
n=1;
y0=0;%循环条件预赋值
y1=0.5*h*（f（a）+f（b））;%循环前的一步
whileabs（y1-y0）>error
sum=0;
xs=a+0.5*h;
y0=y1;%后推
fori=1:
n
sum=sum+f（xs）;
xs=xs+h;
end
y1=0.5*y1+0.5*h*sum;%1式
n=2*n;%保证循环变量才置后，且不引起不一致性
h=h/2;
%y1=y0;%如果与1式合并，无法控制循环
end
I=y1;
End
计算上三式得值：
[I,n]=RecursionTrapezoid（'x^0.5*log（x）',0,1）
I=
NaN
n=
1
都无法求解，前两个是log（0），sin（0）/0无法计算，后者是无穷积分。
采用折中的方法是，二者都在积分区间一致有界，误差可估计，故用1.0e-6代替0求近似解
运行都超过十分钟，终止计算。
采用内联函数inline方法加快效率（代替sym）
第二式I=0.946057561038914n=32I=0.946082070343826n=32768（为1.0e-10）
I=-0.444395950087616n=1024第一式n为插值点数
复合辛普森公式
function[I,step]=IntSimpson（f,a,b,type,eps）
%辛普森系列公式求函数f在区间[a,b]上的定积分
if（type==3&&nargin==4）
disp（'lackvarargin'）;
end
I=0;
switchtype
case1,%%辛普森公式
I=（（b-a）/6）*（f（a）+...
4*f（（a+b）/2）+...
f（b））;
step=1;
case2,%%辛普森3/8公式
I=（（b-a）/8）*（f（a）+...
3*f（（2*a+b）/3）+...
3*f（（a+2*b）/3）+f（b））;
step=1;
case3,%%复合辛普森公式
n=2;
h=（b-a）/2;
I1=0;
I2=（f（a）+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/6）*（f（x）+...
4*f（（x+x1）/2）+...
f（x1））;
end
end
I=I2;
step=n;
end
>>[I,step]=IntSimpson（f,a,b,type,eps）计算前二式得：
I=-0.444332482271752n=164（n为2分次数）I=0.946082310887237n=4
龙贝格公式
function[I,step]=Romberg（f,a,b,error）
ifnargin<3
disp（'Youmustinputthreevariablesatleast'）;
return;
elseifnargin==3
error=1.0e-4;%1.0e-10不现实
end
%h=b-a;
%preT
（1）=0.5*h*（f（a）+f（b））;%可扩增向量
%exitflag=0;
%k=2;
%n=1;%增加的节点数
%while（~exitflag）%该算法则指用到上两行：
对角行和次对角行
%xs=a+0.5*h;
%sum=0.;
%fori=1:
n
%sum=sum+f（xs）;
%xs=xs+h;
%end
%postT
（1）=0.5*preT
（1）+0.5*h*sum;
%fori=2:
k
%postT（i）=4^（i-1）/（4^（i-1）-1）*postT（i-1）-1/（4^（i-1）-1）*preT（i-1）;
%end
%if（abs（postT（k）-preT（k-1））%exitflag=1;
%end
%fori=1:
k
%preT（i）=postT（i）;%前推
%end
%postT（k）
%k=k+1;
%n=2*n;
%h=h/2;
%end
M=1;
tol=10;%误差
k=0;
T=zeros（1,1）;%行表示二分次数，列表示加速次数，只有第一列二分产生，其余均递推。
h=b-a;
T（1,1）=（h/2）*（f（a）+f（b））;%初始值,T可扩增。
whiletol>error
k=k+1;
h=h/2;
Q=0;%sum
fori=1:
M
x=a+h*（2*i-1）;
Q=Q+f（x）;
end
T（k+1,1）=T（k,1）/2+h*Q;%梯形步长减半的复合。
这是第二步
M=2*M;%节点数翻倍
forj=1:
k%这是把整个下三角都求出来。
T（k+1,j+1）=T（k+1,j）+（T（k+1,j）-T（k,j））/（4^j-1）;%简单的变换，分离个1出
%来，输入更方便，下标从1开始故j+1
end
tol=abs（T（k+1,j+1）-T（k,j））;
end
I=T（k+1,k+1）;
step=k;%第一列k次二分。
End
I=-0.444444214475247step=14
I=0.946082070387222step=3
高斯拉盖尔求积公式
functionI=IntGaussLager（f,n,AK,XK）
%高斯拉盖尔，n大于5自加精度,AK,XK,可自定义，默认最多是五项。
if（n<6&&nargin==2）
AK=0;
XK=0;
else
I=sum（AK.*f（XK））;
end
switchn
case2,
I=0.853553*f（-0.585786）+...
0.146447*f（3.414214）;
case3,
I=0.711093*f（0.415575）+...
0.278518*f（2.294280）+...
0.0103893*f（6.289945）;
case4
I=0.603154*f（0.322548）+...
0.357419*f（1.745761）+...
0.0388879*f（4.536620）+...
0.000539295*f（9.395071）;
case5
I=0.521756*f（0.263560）+...
0.398667*f（1.413403）+...
0.0759424*f（3.596426）+...
0.00361176*f（7.085810）+...
0.0000233700*f（12.640801）;
End
求解第三式：
ans=0.498903451386692
高斯列让德公式
functionI=IntGauss（f,a,b,n,AK,XK）
%%AK£¬XKΪ×Ô¶¨ÒåµÄϵÊýÓë×ø±ê¡£´úÊý¾«¶È¿É¼Ó1¡£¸ß˹ÁÐÈõÂ
if（n<6&&nargin==4）
AK=0;
XK=0;
else
XK1=（（b-a）/2）*XK+（（a+b）/2）;
I=（（b-a）/2）*sum（AK.*f（XK1））;
end
ta=（b-a）/2;
tb=（a+b）/2;
switchn
case0,%%...·Ö¿ªf（x0）+...f（x1）
I=2*ta*f（tb）;
case1,
I=ta*（subs（sym（f）,findsym（sym（y））,ta*0.5773503+tb）+...
f（-ta*0.5773503+tb））;
case2,
I=ta*（0.55555556*f（ta*0.7745967+tb）+...
0.55555556*f（-ta*0.7745967+tb）+...
0.88888889*f（tb））;
case3,
I=ta*（0.3478548*f（ta*0.8611363+tb）+...
0.3478548*f（-ta*0.8611363+tb）+...
0.6521452*f（ta*0.3398810+tb）+...
0.6521452*f（-ta*0.3398810+tb））;
case4,
I=ta*（0.2369269*f（ta*0.9061793+tb）+...
0.2369269*f（-ta*0.9061793+tb）+...
0.4786287*f（ta*0.5384693+tb）+...
0.4786287*f（-ta*0.5384693+tb）+...
0.5688889*f（tb））;
case5,
I=ta*（0.1713245*f（ta*0.9324695+tb）+...%%¾ßÓжԳÆÐÔx（k）
0.1713245*f（-ta*0.9324695+tb）+...
0.3607616*f（ta*0.6612094+tb）+...
0.3607616*f（-ta*0.6612094+tb）+...
0.4679139*f（ta*0.2386192+tb）+...
0.4679139*f（-ta*0.2386192+tb））;
end
第一式ans=-0.446183493284934
第二式ans=0.946083069543033
>>loadAg.dat
>>Ag
Ag=
10.000000000000000-7.00000000000000001.000000000000000
-3.0000000000000002.0999990000000006.0000000000000002.000000000000000
5.000000000000000-1.0000000000000005.000000000000000-1.000000000000000
2.0000000000000001.00000000000000002.000000000000000
>>b1=[85.90000151]'
b1=
8.000000000000000
5.900001000000000
5.000000000000000
1.000000000000000
>>GaussOrder（Ag,b1）
ans=
1.0e+007*
0.000000800000000
0.000000830000100
2.075000349709961
0.000000100000000
>>Gauss_pivot（Ag,b1）
ans=
0.000000000000000
-1.000000000000000
1.000000000000000
1.000000000000000
差异：
显然用高斯顺序消去法得到的结果是错误的，用高斯列主元消去法得到了正确解。
造成这种区别的原因可能是：
后者避免了微小量作为列主元（即分母），不会放大误差，保证了算法的稳定性，避免了误差危害
高斯顺序消元法
functionx=GaussOrder（A,b）
ifnargin==2
A=[A,b];
end
N=size（A）;
s=N
（1）;
ifN
（1）~=N
（2）-1&&length（N）~=2
disp（'error-input'）;
return;
end
fori=1:
s%A约化为上三角形矩阵
ifA（i,i）==0
disp（'对角元素有0'）;
return;
end
forj=i+1:
s
A（j,i）=-A（j,i）/A（i,i）;
fork=i+1:
N
（2）
A（j,k）=A（j,k）+A（j,i）*A（i,k）;
end
end
end
functionX=BackSubstitution（A）%%这里用U，b不方便
N=size（A）;
s=N
（1）;
X=ones（s,1）;
A（s,s+1）=A（s,s+1）/A（s,s）;
fori=s-1:
1
forj=s:
i+1
A（i,s+1）=A（i,s+1）-A（j,s+1）*A（i,j）;
end
A（i,s+1）=A（i,s+1）/A（i,i）;
end
fori=1:
s
X（i）=A（i,s+1）;
end
end
x=BackSubstitution（A）;
end
列主元消去法
function[X,det]=Gauss_pivot（A,b）
ifnargin==0
disp（'没参数'）
return;
end
N=size（A）;
n=length（b）;
ifN
（1）~=N
（2）&&length（N）~=2
disp（'输入A有误'）;
return;
end
ifN
（1）~=n
disp（'A和b不匹配'）;
return;
end
det=1;
max=0;
t=0;
c=0;
c2=0;
fori=1:
n
%%选列主元
forj=i:
n%%找i列主元
ifmaxmax=A（j,i）;
t=j;
end
end
ifmax<1.0e-6%%max=0
disp（'没找到列主元'）;
return;
end
%%交换列主元至上行
ift~=i
fork=i:
n
c=A（t,k）;
A（t,k）=A（i,k）;
A（i,k）=c;
end
c2=b（i）;
b（i）=b（t）;
b（t）=c2;
det=-det;
end
%A约化为上三角形矩阵,消元第i列的过程；
ifA（i,i）==0
disp（'对角元素有0，无法继续'）;
return;
end
det=A（i,i）*det;
forj=i+1:
n
A（j,i）=-A（j,i）/A（i,i）;
fork=i+1:
n
A（j,k）=A（j,k）+A（j,i）*A（i,k）;
end
b（j）=b（j）+A（j,i）*b（i）;
end
end
functionX=BackSubstitution（A,b）
N=size（A）;
s=N
（1）;
b（s）=b（s）/A（s,s）;
fori=s-1:
1
forj=s:
i+1
b（i）=b（i）-b（j）*A（i,j）;
end
b（i）=b（i）/A（i,i）;
end
X=b;
end
X=BackSubstitution（A,b）;
end
追赶法
functionx=followup（A,b）
n=rank（A）;
for（i=1:
n）
if（A（i,i）==0）
disp（'Error:
对角元素有0!
'）;
return;
end
end
d=ones（n,1）;
a=ones（n-1,1）;
c=ones（n-1,1）;
for（i=1:
n-1）
a（i,1）=A（i+1,i）;
c（i,1）=A（i,i+1）;
d（i,1）=A（i,i）;
end
d（n,1）=A（n,n）;
%求解Ly=b的解y，解保存在b中
for（i=2:
n）
d（i,1）=d（i,1）-（a（i-1,1）/d（i-1,1））*c（i-1,1）;%%对角d的更新
b（i,1）=b（i,1）-（a（i-1,1）/d（i-1,1））*b（i-1,1）;%%可以存在b里（存在a也可以）。
end
x（n,1）=b（n,1）/d（n,1）;
for（i=（n-1）:
-1:
1）
x（i,1）=（b（i,1）-c（i,1）*x（i+1））/d（i,1）;
end
end
2.用列主元消去法解如下线性方程组
比较两方程组解的差异，并指出产生差异的原因.
>>loadAo1.dat
>>b2=[111]'
b2=
1
1
1
>>Gauss_pivot（Ao1,b2）
ans=
1.0e+003*
1.592599624841381
-0.631911376202549
-0.493617724759390
>>loadAo2.dat
>>Gauss_pivot（Ao2,b2）
ans=
1.0e+002*
1.195273381259593
-0.471426044312964
-0.368402561091259
差异：
相差一个数量级，可见微小的系数差别可能造成结果的巨大差异，这是因为消元过程中不同量得权值在变化，且变化巨大，每次消元都有误差积累。
3.编写追赶法的程序代码，并求解教材第177页的第9题。
>>loadflp.dat
>>followup（flg,b3）
ans=
0.833333333333333
0.666666666666667
0.500000000000000
0.333333333333333
0.166********6667
雅可比迭代法
functiony=Jacobi（A,b）
[m,n]=size（A）;
ifm~=n
disp（'error'）;
quit;
end
ifm~=length（b）
disp（'error'）;
quit;
end
x0=zeros（n,1）;
flag=1;
whileflag
fori=1:
n
x1（i）=（b（i）-A（i,:
）*x0）/A（i,i）+x0（i）;
end
ifnorm（x1.-x0.）<1.0e-7
flag=0;
end
x1=x0;
end
y=x1;
end
高斯塞德尔
functionGaussSeidel（A,b,x0,eps,M）
ifnargin==3
eps=1.0e-6;
M=200;
elseifnargin==4
M=200;
elseifnargin<3
disp（'lackinput-error'）;
return;
end
D=diag（diag（A））;%求A的对角矩阵
L=-tril（A,-1）;
U=-triu（A,1）;
G=（D-L）\U;
f=（D-L）\b;
x=G*x0+f;
n=1;
whilenorm（x-x0）>=eps
x0=x;
x=G*x0+f;
n=n+1;
if（n>=M）
disp（'Warning:
迭代次数太多，可能不收敛'）;
return;
end
end
SOR迭代法
function[x,n]=SOR（A,b,x0,w,eps,M）
ifnargin==4
eps=1.0e-6;
M=200;
elseifnargin<4
disp（'error'）;
return;
elseifnargin==5
M=200;
end
if（w<=0||w>=2）
disp（'error'）;
return;
end
D=diag（diag（A））;
L=-tril（A,-1）;
U=-triu（A,1）;
B=inv（D-L*w）*（（1-w）*D+w*U）;
f=w*inv（（D-L*w））*b;
x=B*x0+f;
n=1;
whilenorm（x-x0）>=eps
x0=x;
x=B*x0+f;
n=n+1;
if（n>=M）
disp（'warning:
迭代次数太多，可能不收敛'）;
return;
end
end
1.编写雅克比迭代法、高斯-赛德尔迭代法的程序代码；取相同的初始点，分别用所得的程序求解如下方程组
比较两种方法所得的结果；如果有差异，解释结果差异的原因。
ans=（初值（0,-1,2）雅可比迭代的两个结果
-0.135********5135ans=-331（初值（000）
.0810********
3.918918918918919
高斯-塞德尔：
1.x=2.x=
-0.135********5135-3
.0810********3
3.9189189189189191
2.用SOR法的程序解如下线性方程组，松弛因子分别为
要求当
时迭代终止。
比较不同采用松弛因子时的迭代次数。
ans=（w取0.9且初值是0向量）ans=（w取1.0且初值是0向量）
-4.000000039257963-4.000000231430727
3.0000000019232012.999999856370402
2.00
升级会员 