matlab程序集.docx
- 文档编号:3665098
- 上传时间:2022-11-24
- 格式:DOCX
- 页数:22
- 大小:20.77KB
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给定一个误差初值以进入循环计算
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlab 程序