1、matlab程序集力学中的保守力场与非保守力场syms x y z thetax=cos(theta);y=sin(theta);z=7*theta;Fx+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(保守立场,做工与路径无关)syms x y z thetax=cos(theta);
2、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)复变函数与积分变换复数与复矩阵的形成a
3、1=7+8*ia2=5*exp(6*i)a3=2+2*i 4-4*i 5+6*i3-5*i 2-2*i 4-8*ib1=randn(4,4);b2=rand(4,4);format shorta4=b1+b2*i复数的基本运算disp(产生正态随即矩阵)a1=randn(4,4)disp(产生hilbert矩阵)a2=hilb(4)disp(生成的a矩阵为)A=a1+a2*idisp(实部为)re=real(A)disp(虚部为)im=imag(A)disp(A的各元素模为)rou=abs(A)disp(A的个元素辐角为)theta=angle(A)disp(A的共轭为)B=conj(A)di
4、sp(A及其共轭的和为)C=A+Bb=3+4*i;disp(A中各元素与b的乘积为)D=b.*Adisp(A中各元素除以b)E=A/bdisp(A中各元素的sin函数为)s=sin(A)disp(A中各元素的三次米)F=A.3disp(A中各元素的对数函数为)H=log(A)disp(A中各元素的指数函数)G=exp(A)注意在matlab中单引号不识别留数s2=5 3 -2 7;s1=-4 0 8 3;r,p,k=residue(s2,s1)syms z a bf=(exp(i*a*z)-exp(i*b*z)/z2c_f=limit(f*z,z,0)g=1/(z-1)2*exp(z2)c_g
5、=limit(diff(g*(z-1)2,z,1)/prod(1:1),z,1)留数在计算曲线积分中的应用syms zf=1/(z2-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)/z3;c3=limit(diff(g*z3),z,2)/prod(1:2),z,0);s2=2*pi*i*c3傅立叶变换syms xf=exp(-2*abs(x)F=fourier(f)subplot(1,2,1)ezplot(f)gridsubplot(1,2,2
6、)ezplot(F)gridsyms wF=(dirac(w-2)+dirac(w+2)/2f=ifourier(F)subplot(1,2,1)ezplot(F)gridsubplot(1,2,2)ezplot(f)grid拉普拉斯变换syms tf=sin(2*t)+sinh(3*t)Lf=laplace(f)pretty(Lf)subplot(1,2,1)ezplot(f)gridsubplot(1,2,2)ezplot(Lf)grid拉普拉斯逆变换syms sLf=1/(s+2)/(s+3)f=ilaplace(Lf)pretty(f)subplot(2,2,1)ezplot(Lf)g
7、ridsubplot(2,2,2)ezplot(f)gridjacobi迭代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=Db;x=B*x0+f;k=1;while norm(x-x0)=tol x0=x; x=B*x0+f; k=k+1; if(k=max1) return; end k xenda=4 -1 1;4 -8 1;-2 1 5;b=7 -21 15;x0=0 0 0;x,k=Fjacobi(a,b,x0,1e-7)Gauss-Seidel
8、迭代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;while norm(x-x0)=tol x0=x; x=G*x0+f; k=k+1; if(k=max1) return; endenda=4 -1 1;4 -8 1;-2 1 5;b=7 -21 15;x0=0 0 0;x,k=Fgseid(a,b,x0,1e-7)逐次超松弛迭代法function x,k=Fsor(A,b,x0,w,tol)max=300;if(
9、w=2)error;return;endD=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;while norm(x-x0)=tol x0=x; x=B*x0+f; k=k+1; if(k=max) return; endenda=5 -1 -1 -1-1 10 -1 -1-1 -1 5 -1-1 -1 -1 10;b=-4 12 8 34;x0=1 1 1 1;x,k=Fsor(a,b,x0,1.2,1e-7)高斯消元法function X=Fgauss
10、(A,b)zg=A b;n=length(b);ra=rank(A);rz=rank(zg);temp1=rz-ra;if temp10,returnendif ra=rz if ra=nX=zeros(n,1);C=zeros(1,n+1);for p= 1:n-1for k=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);endendb=zg(1:n,n+1);A=zg(1:n,1:n);X(n) =b(n)/A(n,n);for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1
11、:n)/A(q,q); end elseendenda=2 2 34 7 7-2 4 5;b=3 1 -7;x=Fgauss(a,b)列主元消去法从这里开使打印function X=Fzhuyuan(A,b)zg=A b;n=length(b);ra=rank(A);rz=rank(zg);temp1=rz-ra;if temp10,disp(方程组无一般意义下的解,系数矩阵与增广矩阵秩不同)returnendif ra=rz if ra=nX=zeros(n,1);C=zeros(1,n+1);for p= 1:n-1Y,j=max(abs(zg(p:n,p);C=zg(p,:);zg(p,
12、:)=zg(j+p-1,:);zg(j+p-1,:)=C;for k=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);endendb=zg(1:n,n+1);A=zg(1:n,1:n);X(n) =b(n)/A(n,n);for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q); end else disp(方程组为欠定方程组)endenda=0 2 0 12 2 3 24 -3 0 16 1 -6 -5;b=0 -2 -7 6;x=Fzhuyuan(a,b)LU
13、分解法计算线性方程组a=5.8 -1 -1 4.67 -8 1 -30.39.5 2 5 -16 -1 12.9 10;b=21.3 -15.7 16.6 7.9;l,u=lu(a);x=u(lb) 注意是字母l,不是数字1CHOLESKY分解法计算线性方程组楚列斯基分解AX=b(系数矩阵)a=16.42904628903813 16.12007848664000 7.82278120608096 2.9570087900000016.12007848664000 17.16712678518433 4.36553354006000 4.363105971000007.822781206080
14、96 4.36553354006000 14.88811554662400 -0.532619408000002.95700879000000 4.36310597100000 -0.53261940800000 18.86911452484754;b=15.6 5.8 20.8 11.9;ch=chol(a);x=ch(chb)奇异值分解法计算线性方程组AX=ba=6.5 -1 -1 3.66.2 7 -5 43 2.1 -6 4.81 5.6 3.7 2.1;b=12.3 21.4 -7.8 21;u,s,v=svd(a)x=v*inv(s)*u*b双共轭梯度法计算AX=bA=23.6 -
15、6.8 4.5 6.3 7.2 8.64.2 7.62 0.8 -3.1 -4.5 5.3 9.8 -6.8 45.2 3.07 -4.37 7.828.4 -9.05 24.4 2.65 6.53 -3.64-9.5 4.54 15.8 4.89 4.8 -7.4826.32 3.08 4.6 -3.5 -7.8 25.6;b=-9.5 8.12 42.3 -21.4 4.7 3.47;%调用函数计算x,flag,relres,iter,resvec=bicg(A,b,1e-12)%x输出运算结果%flag给出计算是否失败信息%relres计算误差范数%实际iter迭代计算次数、%每次迭代误
16、差相对范数plot(resvec)%输出每次迭代误差曲线共轭梯度的LSQR方法计算AX=bA=2.6 -3.1 0.5 0.3 1.2 8.8 3.2 7.6 0.8 -9.1 -4.5 5.3 9.8 -6.4 4.2 3.37 -4.3 3.3 6.4 -4.05 21.4 2.15 6.5 -2.4 -9.5 4.54 15.8 4.89 4.8 -7.4 26.3 3.4 4.5 -3.51 -7.18 2.6 4.6 3.4 4.2 -2.15 10.4 0.22 -1.5 8.2 4.1 -0.41 1.71 2.4 -9.15 4.5 15.8 4.8 4.5 -7.8 8.4
17、-9.05 24.4 2.65 6.53 -3.64 0.2 7.61 10.8 -5.1 -4.6 4.3 18 -0.8 4.2 3.7 -5.7 7.3;b=-6.87.624.543.082.654.89-3.5-2.5-21.4-9.053.63-1.8;%调用函数计算x,flag,relres,iter,resvec=lsqr(A,b,20)%x输出运算结果%flag给出计算是否失败信息%relres计算误差范数%iter迭代次数%resvec每次迭代误差相对范数%默认期望误差为1e-6%允许迭代20次plot(resvec)%作图画出误差范数变化线性方程组的最小残差法(minre
18、s)计算AX=btemp=0.6 -6.8 4.5 6.3 7.2 8.61.2 7.62 0.8 -3.1 -4.5 5.31.8 -6.8 0.2 3.07 -4.37 7.821.4 -9.05 0.4 2.65 0.53 -3.64-9.5 -4.54 -15.8 4.89 4.8 -7.482.32 -3.08 4.6 -3.5 -7.8 05.6;A=temp*temp%产生对称的方阵,系数矩阵A为题目所要求的矩阵b=-6.5 7.8 2.3 -1.4 7.5 4.2;x,flag,relres,iter,resvec=minres(A,b,20)%x返回运算结果%flag判断运算
19、是否成功%relres计算误差范数%iter迭代次数%resvec每次迭代误差相对范数 plot(resvec)%作图画出误差范数变化线性方程组的广义最小残差法(minres)计算AX=bA=25.4 -6.5 5.5 6.3 6.2 4.64.5 7.62 0.7 -3.1 -4.5 5.18.9 -6.8 5.2 3.07 -4.37 7.88.4 -5.05 4.4 2.65 6.53 -4.4-9.5 4.54 6.3 4.89 4.8 -7.46.3 4.08 1.6 -3.5 -7.8 5.6;b=-4.5 7.2 4.3 -2.4 24.7 5.7x,flag,relres,it
20、er,resvec=gmres(A,b)plot(resvec)第五章 非线性方程的求根波尔查诺二分法(区间半分法)使用二分法计算方程f(x)=(x-1)3-3x+2在区间【2,4】上的根,并把通过数值方法绘制函数比较计算结果是否zhengquefunction gen=erfen(f,a,b,tol)%f为方程f(x)=0中的f(x)%如果输入变量缺省,则默认误差精度为1E-3if(nargin=3)tol=1.0e-3;endgen=compute_gen(f,a,b,tol)function r=compute_gen(f,a,b,tol)%计算左端点函数值fa=subs(f,a);%右
21、端函数值fb=subs(f,b);区间中的函数值fzd=subs(f,(a+b)/2);sub函数 R=sub(S,old,new)其中S为符号表达式,old为老变量,new为新变量Sub把老变量替换为新变量if(fa*fzd0) 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)x1=subs(f,x0)+x0;迭代计算 wucha=abs(x1-x0); x0=x1;更新x0的值 在循环中 这一句非常重要 time=time+1;记下迭代次数endx=x1;以p
22、icard.m保存不动点迭代测试主函数x,time=Picard(35*x4-30*x2+3)/8,0.3)x=0:0.01:1;f=(35*x.4-30*x.2+3)/8;plot(x,f)gridtitle(四阶勒让德多项式)Aitken(艾特肯)加速方法计算3阶勒让德多项式=0在增量=0.1附近的根,其中3阶的勒让德多项式为P3(x)=(5x3-3x)/2并做出图形function gen,time=Aitken(func,x0,tol)缺省的情况下,误差容限为是的-5次方if(nargin=2)tol=1.0e-5;endgen=x0;x(1:2)=0,0t=0;记录迭代次数m=0;x
23、2=x0;wucha=0.1;设置误差初值while(wuchatol)t=t+1;记下积累一次迭代次数x1=x2;temp=gen;gen=subs(func,temp)+temp;x(t)=gen;if(t2) 迭代超过两次使用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);endendgen=x2;time=t;记录迭代次数以文件名Aitken.m命名x_Aitken,time_Aitken=Aitken(1/2*(5*x3-3*x),0.1)x_picard,tim
24、e_picard=Picard(1/2*(5*x3-3*x),0.1)x=-1:0.01:1;f=1/2*(5*x.3-3*x);plot(x,f)gridSteffense迭代法非线性方程x2-5*cos (2*x)-4在x=1附近的根。画出函数图像,并根据其与x轴的交点的分布情况对根结果进行分析function gen,time=Steff(fun,x0,tol)缺省的情况下,误差容限为是的-5次方if(nargin=2)tol=1.0e-5;endtime=0;记录迭代次数wucha=0.1设置前后两次迭代的误差gen=x0;while(wuchatol)x1=gen;y=subs(fu
25、n,x1)+x1;z=subs(fun,y)+y;gen=x1-(y-x1)2/(z-2*y+x1);加速公式 wucha=abs(gen-x1); time=time+1;迭代加一次的记录endgen;计算结果time;迭代次数function Steff_main()x_steff,time_steff=Steff(x2-5*cos(2*x)-4,1)x=-4:0.02:4;y=x.2-5*cos(2*x)-4;plot(x,y)gridnewton-raphson(牛顿拉夫森)迭代方法计算方程f(x)=x3+2x2+10x-20在【1,2】区间内的一个根并做出图形分析结果function
26、 gen,time=Newton(f,x0,tol)x0为迭代初值tol为指定误差容限 如果缺省默认为10的-5次方if(nargin=2)tol=1.0e-5;enddf=diff(sym(f);计算原函数的导数x1=x0;time=0;wucha=0.1给定一个误差初值以进入循环计算while(wuchatol) time=time+1;fx=subs(f,x1);df=subs(df,x1);gen=x1-fx/df;迭代公式 wucha=abs(gen-x1); x1=gen;把迭代后的值赋给x1endgen以newton.m命名disp(使用newton-raphson迭代情况)x,
27、time=Newton(x3+2*x2+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;enddf=diff(sym(f);计算原函数的导数df2=diff(df);二阶导数x1=x0;time=0;wucha=0.1给定一个误差初值以进入循环计算
copyright@ 2008-2022 冰豆网网站版权所有
经营许可证编号:鄂ICP备2022015515号-1