MATLAB源代码.docx
《MATLAB源代码.docx》由会员分享，可在线阅读，更多相关《MATLAB源代码.docx（37页珍藏版）》请在冰豆网上搜索。
MATLAB源代码
埃特金插值
functionf=Atken（x,y,x0）
symst;
if（length（x）==length（y））
n=length（x）;
else
disp;
return;
end
y1（1:
n）=t;
for（i=1:
n-1）
for（j=i+1:
n）
y1（j）=y（j）*（t-x（i））/（x（j）-x（i））+y（i）*（t-x（j））/（x（i）-x（j））;
end
y=y1
simplify（y1）;
end
if（nargin==3）
f=subs（y1（n）,'t',x0）;
else
simplify（y1（n））;
f=collect（y1（n））;
f=vpa（f,6）;
end
特征多项式法
functionl=Chapoly（A）
%特征多项式法求矩阵特征值
%已知矩阵：
A
%求得的矩阵特征值：
l
symst;
N=size（A）;
n=N（1,1）;
y=det（A-t*eye（n,n））;
l=solve（y）;
l=vpa（1,5）;%
functionf=Chebyshev（y,k,x0）
symst;
T（1:
k+1）=t;
T
（1）=1;
T
（2）=t;
c（1:
k+1）=0.0;
c
（1）=int（subs（y,findsym（sym（y））,sym（'t'））*T
（1）/sqrt（1-t^2）,t,-1,1）/pi;
c
（2）=2*int（subs（y,findsym（sym（y））,sym（'t'））*T
（2）/sqrt（1-t^2）,t,-1,1）/pi;
f=c
（1）+c
（2）*t;
fori=3:
k+1
T（i）=2*t*T（i-1）-T（i-2）;
c（i）=2*int（subs（y,findsym（sym（y））,sym（'t'））*T（i）/sqrt（1-t^2）,t,-1,1）/2;
f=f+c（i）*T（i）;
f=vpa（f,6）;
if（i==k+1）
if（nargin==3）
f=subs（f,'t',x0）;
else
f=vpa（f,6）;
end
end
end
弦割法求解非线性方程
function[xkt]=ChordsecantToEquation（f,x0,x1,eps）
%弦割法求解非线性方程
%[xkt]=ChordsecantToEquation（f,x0,x1,eps）
%x:
近似解
%k:
迭代次数
%t:
运算时间
%f:
原函数，定义为内联函数
%x0,x1:
初始值
%eps:
误差限
%
%应用举例：
%f=inline（'x^3+4*x^2-10'）;
%x=ChordsecantToEquation（f,1,2,0.5e-6）
%[xk]=ChordsecantToEquation（f,1,2,0.5e-6）
%[xkt]=ChordsecantToEquation（f,1,2,0.5e-6）
%函数的最后一个参数也可以不写，默认情况下，eps=0.5e-6
%[xkt]=ChordsecantToEquation（f,1,2）
ifnargin==3
eps=0.5e-6;
end
tic;
k=0;
while1
x=x1-f（x1）*（x1-x0）./（f（x1）-f（x0））;
k=k+1;
ifabs（x-x1）30
break;
end
x0=x1;
x1=x;
end
t=toc;
ifk>=30
disp（'迭代次数太多。
'）;
x=0;
t=0;
end
复合求积公式
function[I,step]=CombineTraprl（f,a,b,eps）
%f被积函数
%a,b积分上下限
%eps精度
%I积分结果
%step积分的子区间数
if（nargin==3）
eps=1.0e-4;
end
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
I=I2;
step=n;
functionf=DCS（x,y,x0）
symst;
if（length（x）==length（y））
n=length（x）;
c（1:
n）=0.0;
else
disp（'x和y的维数不相等！
'）;
return;
end
c
（1）=y
（1）;
for（i=1:
n-1）
for（j=i+1:
n）
y1（j）=（x（j）-x（i））/（y（j）-y（i））;
end
c（i+1）=y1（i+1）;
y=y1;
end
f=c（n）;
for（i=1:
n-1）
f=c（n-i）+（t-x（n-i））/f;
f=vpa（f,6）;
if（i==n-1）
if（nargin==3）
f=subs（f,'t',x0）;
else
f=vpa（f,6）;
end
end
end;
用欧拉法求一阶常微分方程的数值解
functiony=DEEuler（f,h,a,b,y0,varvec）
formatlong;
N=（b-a）/h;
y=zeros（N+1,1）;
x=a:
h:
b;
y
（1）=y0;
fori=2:
N+1
y（i）=y（i-1）+h*Funval（f,varvec,[x（i-1）,y（i-1）]）;
end
formatshort;
用隐式欧拉法求一阶常微分方程的数值解
functiony=DEimpEuler（f,h,a,b,y0,varvec）
formatlong;
N=（b-a）/h;
y=zeros（N+1,1）;
y
（1）=y0;
x=a:
h:
b;
var=findsym（f）;
fori=2:
N+1
fx=Funval（f,var
（1）,x（i））;
gx=y（i-1）+h*fx-varvec
（2）;
y（i）=NewtonRoot（gx,-10,10,eps）;
end
formatshort;
龙格库塔法
functiony=DELGKT2_suen（f,h,a,b,y0,varvec）
formatlong;
N=uint16（（b-a）/h）;
y=zeros（N+1,1）;
y
（1）=y0;
x=a:
h:
b;
var=findsym（f）;
fori=2:
N+1
K1=Funval（f,varvec,[x（i-1）y（i-1）]）;
K2=Funval（f,varvec,[x（i-1）+2*h/3y（i-1）+K1*2*h/3]）;
y（i）=y（i-1）+h*（K1+3*K2）/4;
end
formatshort;...
用改进欧拉法求一阶常微分方程的数值解
functiony=DEModifEuler（f,h,a,b,y0,varvec）
formatlong;
N=（b-a）/h;
y=zeros（N+1,1）;
y
（1）=y0;
x=a:
h:
b;
var=findsym（f）;
fori=2:
N+1
v1=Funval（f,varvec,[x（i-1）y（i-1）]）;
t=y（i-1）+h*v1;
v2=Funval（f,varvec,[x（i）t]）;
y（i）=y（i-1）+h*（v1+v2）/2;
end
formatshort;
一阶微分方程
Functiondy=myfun_1（x,y）%定义输入,输出变量和函数文件名
dy=zeros（2,1）;%明确dy的维数,用微分方程组时不可缺省
dy
（1）=y
（2）;%dy（m）表示y的m阶导数;y（n）表示y的第n列
dy
（2）=2*y
（2）-2*y
（1）+5*exp（2*x）*sin（x）;%与方程组中第二个微分方程对应
基本Guass消去法求解线性方程组
functionx=EqtsBasicGuass（A,b）
%基本Guass消去法求解线性方程组Ax=b
%x=EqtsBasicGuass（A,b）
%x:
解向量，列向量
%A:
线性方程组的矩阵
%b:
列向量
%
%检查输入参数
ifsize（A,1）~=size（b,1）
disp（'输入参数有误！
'）;
x='';
return;
end
%（A|b）
A=[Ab];
%消去过程
n=size（A,1）;
l=zeros（n）;
fork=1:
n-1
fori=k+1:
n
l（i,k）=A（i,k）/A（k,k）;
end
fori=k+1:
n
forj=k+1:
n+1
A（i,j）=A（i,j）-l（i,k）*A（k,j）;
end
forj=1:
k
A（i,j）=0;
end
end
end
%回代过程
x=zeros（n,1）;
x（n）=A（n,n+1）/A（n,n）;
fori=n-1:
-1:
1
y=0;
forj=i+1:
n
y=y+A（i,j）*x（j）;
end
x（i）=（A（i,n+1）-y）/A（i,i）;
end
return;
追赶法求解三对角线性方程组
functionx=EqtsForwardAndBackward（L,D,U,b）
%追赶法求解三对角线性方程组Ax=b
%x=EqtsForwardAndBackward（L,D,U,b）
%x:
三对角线性方程组的解
%L:
三对角矩阵的下对角线，行向量
%D:
三对角矩阵的对角线，行向量
%U:
三对角矩阵的上对角线，行向量
%b:
线性方程组Ax=b中的b，列向量
%
%应用举例:
%L=[-1-2-3];D=[2345];U=[-1-2-3];b=[61-21]';
%x=EqtsForwardAndBackward（L,D,U,b）
%检查参数的输入是否正确
n=length（D）;m=length（b）;
n1=length（L）;n2=length（U）;
ifn-n1~=1||n-n2~=1||n~=m
disp（'输入参数有误！
'）
x='';
return;
end
%追的过程
fori=2:
n
L（i-1）=L（i-1）/D（i-1）;
D（i）=D（i）-L（i-1）*U（i-1）;
end
x=zeros（n,1）;
x
（1）=b
（1）;
fori=2:
n
x（i）=b（i）-L（i-1）*x（i-1）;
end
%赶的过程
x（n）=x（n）/D（n）;
fori=n-1:
-1:
1
x（i）=（x（i）-U（i）*x（i+1））/D（i）;
end
return;
Guass-Seidel迭代法
function[xk]=EqtsGS（A,b,x0,eps）
%Guass-Seidel迭代法求解线性方程组Ax=b
%[xk]=EqtsGS（A,b,x0,eps）
%x:
解向量，列向量
%k:
迭代次数
%A:
系数矩阵
%b:
列向量
%x0:
迭代初始值，列向量
%eps:
误差限，可缺省，缺省值为0.5e-6
%
%应用举例:
%A=[1031;2-103;1310];b=[14;-5;14];x0=[0;0;0];
%[xk]=EqtsGS（A,b,x0,0.5e-6）
%x=EqtsGS（A,b,x0）
ifnargin==3
eps=0.5e-6;
end
%检查输入参数
n=length（b）;
ifsize（A,1）~=n||n~=length（x0）
disp（'输入参数有误！
'）;
x='';
k='';
return;
end
%迭代求解
k=0;
x=zeros（n,1）;
while1
k=k+1;
fori=1:
n
z=0;
forj=1:
i-1
z=z+A（i,j）*x（j）;
end
forj=i+1:
n
z=z+A（i,j）*x0（j）;
end
x（i）=（b（i）-z）/A（i,i）;
end
ifnorm（x-x0）<=eps||k>30
break;
end
x0=x;
end
ifk>=30
disp（'迭代次数太多！
'）
x='';
return;
end
return;
Jacobi迭代法求解线性方程组
function[xk]=EqtsJacobi（A,b,x0,eps）
%Jacobi迭代法求解线性方程组Ax=b
%[xk]=EqtsJacobi（A,b,x0,eps）
%x:
解向量，列向量
%k:
迭代次数
%A:
系数矩阵
%b:
列向量
%x0:
迭代初始值，列向量
%eps:
误差限，可缺省，缺省值为0.5e-6
%
%应用举例:
%A=[1031;2-103;1310];b=[14;-5;14];x0=[0;0;0];
%[xk]=EqtsJacobi（A,b,x0,0.5e-6）
%x=EqtsJacobi（A,b,x0）
ifnargin==3
eps=0.5e-6;
end
%检查输入参数
n=length（b）;
ifsize（A,1）~=n||n~=length（x0）
disp（'输入参数有误！
'）;
x='';
k='';
return;
end
%迭代求解
k=0;
x=zeros（n,1）;
while1
k=k+1;
fori=1:
n
z=0;
forj=1:
i-1
z=z+A（i,j）*x0（j）;
end
forj=i+1:
n
z=z+A（i,j）*x0（j）;
end
x（i）=（b（i）-z）/A（i,i）;
end
ifnorm（x-x0）<=eps||k>30
break;
end
x0=x;
end
ifk>=30
disp（'迭代次数太多！
'）
x='';
return;
end
return;
逐次超松驰迭代法
function[xk]=EqtsSOR（A,b,x0,omiga,eps）
%超松弛（SOR,SuccessiveOver-Relaxation）迭代法求解线性方程组Ax=b
%[xk]=EqtsSOR（A,b,x0,eps）
%x:
解向量，列向量
%k:
迭代次数
%A:
系数矩阵
%b:
列向量
%x0:
迭代初始值，列向量
%omiga:
松弛因子，可缺省，缺省值为1，即为GS迭代法
%eps:
误差限，可缺省，缺省值为0.5e-6
%
%应用举例:
%A=[430;34-1;0-14];b=[24;30;-24];x0=[1;1;1];omiga=1.25;
%[xk]=EqtsSOR（A,b,x0,omiga,0.5e-6）
%x=EqtsSOR（A,b,x0）
ifnargin==4
eps=0.5e-6;
end
ifnargin==3
omiga=1;
eps=0.5e-6;
end
%检查输入参数
n=length（b）;
ifsize（A,1）~=n||n~=length（x0）
disp（'输入参数有误！
'）;
x='';
k='';
return;
end
%迭代求解
k=0;
x=zeros（n,1）;
while1
k=k+1;
fori=1:
n
z=0;
forj=1:
i-1
z=z+A（i,j）*x（j）;
end
forj=i+1:
n
z=z+A（i,j）*x0（j）;
end
x（i）=（1-omiga）*x0（i）+omiga*（b（i）-z）/A（i,i）;
end
ifnorm（x-x0）<=eps||k>30
break;
end
x0=x;
end
ifk>=30
disp（'迭代次数太多！
'）
x='';
return;
end
return;
连分式插值
functionf=DCS（x,y,x0）
symst;
if（length（x）==length（y））
n=length（x）;
c（1:
n）=0.0;
else
disp;
return;
end
c
（1）=y
（1）;
for（i=1:
n-1）
for（j=1+1:
n）
y1（j）=（x（j）-x（i）/（y（j）-y（i））;
end
c（i+1）=y1（i+1）;
y=y1;
end
f=c（n）;
for（i=1:
n-1）
f=c（n-i）+（t-x（n-i））/f;
f=vpa（f,6）;
if（i==n-1）
if（nargin==3）
f=subs（f,'t',x0）;
else
f=vpa（f,6）;
end
end
end;
用二分法求方程的一个根
functionroot=HalfInterval（f,a,b,eps）
if（nargin==3）
eps=1.0e-4;
end
f1=subs（sym（f）,findsym（sym（f））,a）;
f2=subs（sym（f）,findsym（sym（f））,b）;
if（f1==0）
root=a;
end
if（f2==0）
root=b;
end
if（f1*f2>0）
disp（'两端点函数值乘积大于0!
'）;
return;
else
root=FindRoots（f,a,b,eps）;
end
functionr=FindRoots（f,a,b,eps）
f_1=subs（sym（f）,findsym（sym（f））,a）;
f_2=subs（sym（f）,findsym（sym（f））,b）;
mf=subs（sym（f）,findsym（sym（f））,（a+b）/2）;
if（f_1*mf>0）
t=（a+b）/2;
r=FindRoots（f,t,b,eps）;
else
if（f_1*mf==0）
r=（a+b）/2;
else
if（abs（b-a）<=eps）
r=（b+3*a）/4;
else
s=（a+b）/2;
r=FindRoots（f,a,s,eps）;
end
end
end
Hermite插值的matlab代码
functionf=Hermite（x,y,y_1,x0）
symst;
f=0.0;
if（length（x）==length（y））
if（length（y）==length（y_1））
n=length（x）;
else
disp（'y和y的导数的维数不相等！
'）;
return;
end
else
disp（'x和y的维数不相等！
'）;
return;
end
fori=1:
n
h=1.0;
a=0.0;
forj=1:
n
if（j~=i）
h=h*（t-x（j））^2/（（x（i）-x（j））^2）;
a=a+1/（x（i）-x（j））;
end
end
f=f+h*（（x（i）-t）*（2*a*y（i）-y_1（i））+y（i））;
if（i==n）
if（nargin==4）
f=subs（f,'t',x0）;
else
f=vpa（f,6）;
end
end
end
辛普生公式的求积
function[Sn]=IntCompSimpson（f,D4f,a,b,eps）
ifnargin==4
eps=0.5e-7;%默认精度
end
%求被积函数的四次导数在区间[a,b]上的最大值
[x,fval]=fminbnd（D4f,a,b,optimset（'TolX',eps/10））;
fmax=-fval;
%计算等分区间数
n=ceil（sqrt（sqrt（abs（（b-a）^5*fmax/16/180/eps））））;
h=（b-a）/n;%步长
S=f（a）+f（b）+4*f（a+h/2）;
fork=1:
n-1
x1=a+k*h;
x2=x1+h/2;
S=S+4*f（x2）+2*f（x1）;
end
S=S*h/6;
return;
复化梯形公式求积分
function[Tn]=IntCompTrape（f,D2f,a,b,eps）
%复化梯形公式求积分
%[Tn]=IntCompTrape（f,D2f,a,b,eps）
%T:
求积结果
%n:
区间等分数
%f:
被积函数，可利用函数脚本文件事先定义，也可以利用内联函数
?
f:
被积函数的二次导数的相反值，可利用函数脚本文件事先定义，也可以利用内联函数
%取相反值是为了便于计算被积函数的二次导数在区间[a,b]上的最大值
%a:
积分下限
%b:
积分上限
%eps:
误差限
%
%应用举例:
%事先定义
%[Tn]=IntCompTrape（@IntF,@IntD2F,0,1）
%[Tn]=IntCompTrape（@IntF,@IntD2F,0,1,0.5e-7）
%利用内联函数
%F=inline（'3^x'）;D2F=inline（'-（log（3））^2*3^x'）;
%[Tn]=IntCompTrape（F,D2F,0,1）
ifnargin==4
eps=0.5e-7;%默认精度
end
%求被积函数的二次导数在区间[a,b]上的最大值
[x,fval]=fminbnd（D2f,a,b,optimset（'TolX',eps/10））;
fmax=-fval;
%计算等分区间数
n=ceil（sqrt（abs（fmax*（b-a）^3/12/eps）））;
h=（b-a）/n;%步长
T=f（a）+f（b）;
fork=1:
n-1
x1=a+k*h;
T=T+2*f（x1）;
end
T=h*T/2;
return;
拉格朗日插值
（1）
functionf=InterpLagra
升级会员 