matlab程序集.docx
《matlab程序集.docx》由会员分享，可在线阅读，更多相关《matlab程序集.docx（22页珍藏版）》请在冰豆网上搜索。
matlab程序集
力学中的保守力场与非保守力场
symsxyztheta
x=cos（theta）;
y=sin（theta）;
z=7*theta;
F[x+2*y+z+5,2*x+y+z,x+y+z-6];
ds=[diff（x）,diff（y）,diff（z）]';
int_1=int（F*ds,theta,0,2*pi）
x=cos（8*theta）;
y=20*sin（theta）;
z=7*theta;
ds=[diff（x）,diff（y）,diff（z）]';
int_2=int（F*ds,theta,0,2*pi）
disp（'保守立场，做工与路径无关'）
symsxyztheta
x=cos（theta）;
y=sin（theta）;
z=7*theta;
F=[2*x-3*y+4*z-5,z-x+8,x+y+z+12];
ds=[diff（x）,diff（y）,diff（z）]';
int_3=int（F*ds,theta,0,2*pi）
ezplot3（x,y,z,[0,2*pi]）
x=cos（8*theta）;
y=20*sin（theta）;
z=7*theta;
ds=[diff（x）,diff（y）,diff（z）]';
int_4=int（F*ds,theta,0,2*pi）
disp（'非保守立场，做工与路径有关，不存在势能'）
figure
（2）
ezplot3（x,y,z,[0,2*pi]）
复变函数与积分变换
复数与复矩阵的形成
a1=7+8*i
a2=5*exp（6*i）
a3=[2+2*i4-4*i5+6*i
3-5*i2-2*i4-8*i]
b1=randn（4,4）;
b2=rand（4,4）;
formatshort
a4=b1+b2*i
复数的基本运算
disp（‘产生正态随即矩阵’）
a1=randn（4,4）
disp（‘产生hilbert矩阵’）
a2=hilb（4）
disp（‘生成的a矩阵为’）
A=a1+a2*i
disp（‘实部为’）
re=real（A）
disp（‘虚部为’）
im=imag（A）
disp（‘A的各元素模为’）
rou=abs（A）
disp（‘A的个元素辐角为’）
theta=angle（A）
disp（‘A的共轭为’）
B=conj（A）
disp（‘A及其共轭的和为’）
C=A+B
b=3+4*i;
disp（‘A中各元素与b的乘积为’）
D=b.*A
disp（‘A中各元素除以b’）
E=A/b
disp（‘A中各元素的sin函数为’）
s=sin（A）
disp（‘A中各元素的三次米’）
F=A.^3
disp（‘A中各元素的对数函数为’）
H=log（A）
disp（‘A中各元素的指数函数’）
G=exp（A）
注意在matlab中单引号不识别
留数
s2=[53-27];
s1=[-4083];
[r,p,k]=residue（s2,s1）
symszab
f=（exp（i*a*z）-exp（i*b*z））/z^2
c_f=limit（f*z,z,0）
g=1/（z-1）^2*exp（z^2）
c_g=limit（diff（g*（z-1）^2,z,1）/prod（1:
1）,z,1）
留数在计算曲线积分中的应用
symsz
f=1/（z^2-1）*sin（pi*z/4）
c1=limit（f*（z-1）,z,1）;
c2=limit（f*（z+1）,z,-1）;
disp（'封闭曲线s1积分为'）
s1=2*pi*i*（c1+c2）
g=exp（z）/z^3;
c3=limit（diff（（g*z^3）,z,2）/prod（1:
2）,z,0）;
s2=2*pi*i*c3
傅立叶变换
symsx
f=exp（-2*abs（x））
F=fourier（f）
subplot（1,2,1）
ezplot（f）
grid
subplot（1,2,2）
ezplot（F）
grid
symsw
F=（dirac（w-2）+dirac（w+2））/2
f=ifourier（F）
subplot（1,2,1）
ezplot（F）
grid
subplot（1,2,2）
ezplot（f）
grid
拉普拉斯变换
symst
f=sin（2*t）+sinh（3*t）
Lf=laplace（f）
pretty（Lf）
subplot（1,2,1）
ezplot（f）
grid
subplot（1,2,2）
ezplot（Lf）
grid
拉普拉斯逆变换
symss
Lf=1/（s+2）/（s+3）
f=ilaplace（Lf）
pretty（f）
subplot（2,2,1）
ezplot（Lf）
grid
subplot（2,2,2）
ezplot（f）
grid
jacobi迭代
function[x,k]=Fjacobi（A,b,x0,tol）
max1=300;
D=diag（diag（A））;
L=-tril（A,-1）;
U=-triu（A,1）;
B=D\（L+U）;
f=D\b;
x=B*x0+f;
k=1;
whilenorm（x-x0）>=tol
x0=x;
x=B*x0+f;
k=k+1;
if（k>=max1）
return;
end
[kx']
end
a=[4-11;4-81;-215];
b=[7-2115]';
x0=[000]';
[x,k]=Fjacobi（a,b,x0,1e-7）
Gauss-Seidel迭代
function[x,k]=Fgseid（A,b,x0,tol）
max1=300;
D=diag（diag（A））;
L=-tril（A,-1）;
U=-triu（A,1）;
G=（D-L）\U;
f=（D-L）\b;
x=G*x0+f;
k=1;
whilenorm（x-x0）>=tol
x0=x;
x=G*x0+f;
k=k+1;
if（k>=max1）
return;
end
end
a=[4-11;4-81;-215];
b=[7-2115]';
x0=[000]';
[x,k]=Fgseid（a,b,x0,1e-7）
逐次超松弛迭代法
function[x,k]=Fsor（A,b,x0,w,tol）
max=300;
if（w<=0||w>=2）
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））;
x=B*x0+f;
k=1;
whilenorm（x-x0）>=tol
x0=x;
x=B*x0+f;
k=k+1;
if（k>=max）
return;
end
end
a=[5-1-1-1
-110-1-1
-1-15-1
-1-1-110];
b=[-412834]';
x0=[1111]';
[x,k]=Fsor（a,b,x0,1.2,1e-7）
高斯消元法
functionX=Fgauss（A,b）
zg=[Ab];n=length（b）;
ra=rank（A）;
rz=rank（zg）;
temp1=rz-ra;
iftemp1>0,
return
end
ifra==rz
ifra==n
X=zeros（n,1）;
C=zeros（1,n+1）;
forp=1:
n-1
fork=p+1:
n
m=zg（k,p）/zg（p,p）;
zg（k,p:
n+1）=zg（k,p:
n+1）-m*zg（p,p:
n+1）;
end
end
b=zg（1:
n,n+1）;
A=zg（1:
n,1:
n）;
X（n）=b（n）/A（n,n）;
forq=n-1:
-1:
1
X（q）=（b（q）-sum（A（q,q+1:
n）*X（q+1:
n）））/A（q,q）;
end
else
end
end
a=[223
477
-245];
b=[31-7]’;
x=Fgauss（a,b）
列主元消去法从这里开使打印
functionX=Fzhuyuan（A,b）
zg=[Ab];n=length（b）;
ra=rank（A）;
rz=rank（zg）;
temp1=rz-ra;
iftemp1>0,
disp（‘方程组无一般意义下的解，系数矩阵与增广矩阵秩不同’）
return
end
ifra==rz
ifra==n
X=zeros（n,1）;
C=zeros（1,n+1）;
forp=1:
n-1
[Y,j]=max（abs（zg（p:
n,p）））;C=zg（p,:
）;
zg（p,:
）=zg（j+p-1,:
）;zg（j+p-1,:
）=C;
fork=p+1:
n
m=zg（k,p）/zg（p,p）;
zg（k,p:
n+1）=zg（k,p:
n+1）-m*zg（p,p:
n+1）;
end
end
b=zg（1:
n,n+1）;
A=zg（1:
n,1:
n）;
X（n）=b（n）/A（n,n）;
forq=n-1:
-1:
1
X（q）=（b（q）-sum（A（q,q+1:
n）*X（q+1:
n）））/A（q,q）;
end
else
disp（'方程组为欠定方程组'）
end
end
a=[0201
2232
4-301
61-6-5];
b=[0-2-76]’;
x=Fzhuyuan（a,b）
LU分解法计算线性方程组
a=[5.8-1-14.6
7-81-30.3
9.525-1
6-112.910];
b=[21.3-15.716.67.9]’;
[l,u]=lu（a）;
x=u\（l\b）注意是字母l，不是数字1
CHOLESKY分解法计算线性方程组楚列斯基分解AX=b
（系数矩阵）a=[16.4290462890381316.120078486640007.822781206080962.95700879000000
16.1200784866400017.167126785184334.365533540060004.36310597100000
7.822781206080964.3655335400600014.88811554662400-0.53261940800000
2.957008790000004.36310597100000-0.5326194080000018.86911452484754];
b=[15.65.820.811.9]’;
ch=chol（a）;
x=ch\（ch’\b）
奇异值分解法计算线性方程组AX=b
a=[6.5-1-13.6
6.27-54
32.1-64.8
15.63.72.1];
b=[12.321.4-7.821]’;
[u,s,v]=svd（a）
x=v*inv（s）*u’*b
双共轭梯度法计算AX=b
A=[23.6-6.84.56.37.28.6
4.27.620.8-3.1-4.55.3
9.8-6.845.23.07-4.377.82
8.4-9.0524.42.656.53-3.64
-9.54.5415.84.894.8-7.48
26.323.084.6-3.5-7.825.6];
b=[-9.58.1242.3-21.44.73.47]’;
%调用函数计算
[x,flag,relres,iter,resvec]=bicg（A,b,1e-12）
%x输出运算结果
%flag给出计算是否失败信息
%relres计算误差范数
%实际iter迭代计算次数、
%每次迭代误差相对范数
plot（resvec）
%输出每次迭代误差曲线
共轭梯度的LSQR方法计算AX=b
A=[2.6-3.10.50.31.28.8
3.27.60.8-9.1-4.55.3
9.8-6.44.23.37-4.33.3
6.4-4.0521.42.156.5-2.4
-9.54.5415.84.894.8-7.4
26.33.44.5-3.51-7.182.6
4.63.44.2-2.1510.40.22
-1.58.24.1-0.411.712.4
-9.154.515.84.84.5-7.8
8.4-9.0524.42.656.53-3.64
0.27.6110.8-5.1-4.64.3
18-0.84.23.7-5.77.3];
b=[-6.8
7.62
4.54
3.08
2.65
4.89
-3.5
-2.5
-21.4
-9.05
3.63
-1.8];
%调用函数计算
[x,flag,relres,iter,resvec]=lsqr（A,b,[],20）
%x输出运算结果
%flag给出计算是否失败信息
%relres计算误差范数
%iter迭代次数
%resvec每次迭代误差相对范数
%默认期望误差为1e-6
%允许迭代20次
plot（resvec）
%作图画出误差范数变化
线性方程组的最小残差法（minres）计算AX=b
temp=[0.6-6.84.56.37.28.6
1.27.620.8-3.1-4.55.3
1.8-6.80.23.07-4.377.82
1.4-9.050.42.650.53-3.64
-9.5-4.54-15.84.894.8-7.48
2.32-3.084.6-3.5-7.805.6];
A=temp*temp’%产生对称的方阵，系数矩阵A为题目所要求的矩阵
b=[-6.57.82.3-1.47.54.2]’;
[x,flag,relres,iter,resvec]=minres（A,b,[],20）
%x返回运算结果
%flag判断运算是否成功
%relres计算误差范数
%iter迭代次数
%resvec每次迭代误差相对范数
plot（resvec）
%作图画出误差范数变化
线性方程组的广义最小残差法（minres）计算AX=b
A=[25.4-6.55.56.36.24.6
4.57.620.7-3.1-4.55.1
8.9-6.85.23.07-4.377.8
8.4-5.054.42.656.53-4.4
-9.54.546.34.894.8-7.4
6.34.081.6-3.5-7.85.6];
b=[-4.57.24.3-2.424.75.7]’
[x,flag,relres,iter,resvec]=gmres（A,b）
plot（resvec）
第五章非线性方程的求根
波尔查诺二分法（区间半分法）
使用二分法计算方程f（x）=（x-1）^3-3x+2在区间【2，4】上的根，并把通过数值方法绘制函数比较计算结果是否zhengque
functiongen=erfen（f,a,b,tol）
%f为方程f（x）=0中的f（x）
%如果输入变量缺省，则默认误差精度为1E-3
if（nargin==3）
tol=1.0e-3;
end
gen=compute_gen（f,a,b,tol）
functionr=compute_gen（f,a,b,tol）
%计算左端点函数值
fa=subs（f,a）;
%右端函数值
fb=subs（f,b）;
区间中的函数值
fzd=subs（f,（a+b）/2）;
sub函数R=sub（S,old,new）
其中S为符号表达式，old为老变量，new为新变量
Sub把老变量替换为新变量
if（fa*fzd>0）
t=（a+b）/2;
采用递归方法
r=compute_gen（f,t,b,tol）;
else
if（fa*fzd==0）
r=（a+b）/2
else
if（abs（b-a）<=tol）
r=（b+3*a）/4;
else
s=（a+b）/2;
结果
r=compute_gen（f,a,s,tol）;
end
end
end
保存为erfen.m文件
编写检验文件
functionerfen_main（）
x=erfen（'（x-1）^3-3*x+2',2,4）
x=-1:
0.01:
4;
f=（x-1）^3-3*x+2;
plot（x,f）
gridon
保存为erfen_main.m文件
不动点迭代法（Picard迭代）
计算四阶勒让德多项式=0在增量=0.3附近的根，其中四阶的勒让德多项式为
P4（x）=（35x^4-30x^2+3）/8并做出图形
function[x,time]=Picard（f,x0,tol）
结果给出迭代次数
X0为迭代初值
Tol为误差容限
if（nargin==2）
tol=1.0e-5;
end
缺省的情况下，误差容限为是的-5次方
wucha=0.5;设置误差初值
x1=x0;x1=x0为前后两次计算结果
time=0;用于计算迭代次数
while（wucha>tol）
x1=subs（f,x0）+x0;
迭代计算
wucha=abs（x1-x0）;
x0=x1;更新x0的值在循环中这一句非常重要
time=time+1;记下迭代次数
end
x=x1;
以picard.m保存
不动点迭代测试主函数
[x,time]=Picard（'（35*x^4-30*x^2+3）/8',0.3）
x=0:
0.01:
1;
f=（35*x.^4-30*x.^2+3）/8;
plot（x,f）
grid
title（'四阶勒让德多项式'）
Aitken（艾特肯）加速方法
计算3阶勒让德多项式=0在增量=0.1附近的根，其中3阶的勒让德多项式为
P3（x）=（5x^3-3x）/2并做出图形
function[gen,time]=Aitken（func,x0,tol）
缺省的情况下，误差容限为是的-5次方
if（nargin==2）
tol=1.0e-5;
end
gen=x0;
x（1:
2）=[0,0]
t=0;记录迭代次数
m=0;
x2=x0;
wucha=0.1;设置误差初值
while（wucha>tol）
t=t+1;记下积累一次迭代次数
x1=x2;
temp=gen;
gen=subs（func,temp）+temp;
x（t）=gen;
if（t>2）迭代超过两次使用Aitken加速
m=m+1;
x2=x（m）-（x（m+1）-x（m））^2/（x（m+2）-2*x（m+1）+x（m））;
给出两次迭代误差
wucha=abs（x2-x1）;
end
end
gen=x2;
time=t;记录迭代次数
以文件名Aitken.m命名
[x_Aitken,time_Aitken]=Aitken（'1/2*（5*x^3-3*x）',0.1）
[x_picard,time_picard]=Picard（'1/2*（5*x^3-3*x）',0.1）
x=-1:
0.01:
1;
f=1/2*（5*x.^3-3*x）;
plot（x,f）
grid
Steffense迭代法
非线性方程x^2-5*cos（2*x）-4在x=1附近的根。
画出函数图像，并根据其与x轴的交点的分布情况对根结果进行分析
function[gen,time]=Steff（fun,x0,tol）
缺省的情况下，误差容限为是的-5次方
if（nargin==2）
tol=1.0e-5;
end
time=0;记录迭代次数
wucha=0.1设置前后两次迭代的误差
gen=x0;
while（wucha>tol）
x1=gen;
y=subs（fun,x1）+x1;
z=subs（fun,y）+y;
gen=x1-（y-x1）^2/（z-2*y+x1）;加速公式
wucha=abs（gen-x1）;
time=time+1;迭代加一次的记录
end
gen；计算结果
time;迭代次数
functionSteff_main（）
[x_steff,time_steff]=Steff（'x^2-5*cos（2*x）-4',1）
x=-4:
0.02:
4;
y=x.^2-5*cos（2*x）-4;
plot（x,y）
grid
newton-raphson（牛顿拉夫森）迭代方法
计算方程f（x）=x3+2x2+10x-20在【1，2】区间内的一个根并做出图形分析结果
function[gen,time]=Newton（f,x0,tol）
x0为迭代初值
tol为指定误差容限如果缺省默认为10的-5次方
if（nargin==2）
tol=1.0e-5;
end
df=diff（sym（f））;计算原函数的导数
x1=x0;
time=0;
wucha=0.1给定一个误差初值以进入循环计算
while（wucha>tol）
time=time+1;
fx=subs（f,x1）;
df=subs（df,x1）;
gen=x1-fx/df;迭代公式
wucha=abs（gen-x1）;
x1=gen;把迭代后的值赋给x1
end
gen
以newton.m命名
disp（‘使用newton-raphson迭代情况‘）
[x,time]=Newton（'x^3+2*x^2+10*x-20',1.5,1e-4）
画出函数图像与x轴交点的情况
x=0:
0.01:
2;
y=x.^3+2*x.^2+10*x-20;
plot（x,y）
grid
以newton_main.m命名
重根的加速迭代问题
计算方程f（x）=x4-4x2+4=0在x=1.5附近的根，并做出图形分析
function[gen,time]=Newton（f,x0,tol）
x0为迭代初值
tol为指定误差容限如果缺省默认为10的-5次方
if（nargin==2）
tol=1.0e-5;
end
df=diff（sym（f））;计算原函数的导数
df2=diff（df）;二阶导数
x1=x0;
time=0;
wucha=0.1给定一个误差初值以进入循环计算
升级会员 