matlab16常用计算方法

常用计算方法1.超越方程的求解一超越方程为x(2lnx-3)-100=0求超越方程的解。［算法］方法一：用迭代算法。将方程改为100X= 2ln(X0)-3其中x0是一个初始值，由此计算终值x。取最大误差为e=10-4，当|x-x0l>e时，就用x的值换成x0的值，重新进行计算；否则lx-x0l<e为止。［程序］P1_1abs.m如下。%超越方程的迭代算法%清除变量%初始值%空向量%清除变量%初始值%空向量%无限循环%迭代运算%连接结果%如果项数太多则退出循环(暗示发散)%当精度足够高时退出循环%替换初值%结束循环%创建图形窗口xx=[];while1x=100/(2*log(x0)-3);xx=[xx,x];iflength(xx)>1000,break,endifabs(x0-x)<1e-4,break,endx0=x;endfigureplot(xx,'.-','LineWidth',2,'MarkerSize',12)%画迭代线'.-'表示每个点用.来表示，再用线连接gridon %加网格fs=16; %字体大小title('超越方程的迭代折线’,'fontsize',fs)%标题xlabel('\itn', 'fontsize',fs) %x标签ylabel('\itx', 'fontsize',fs) %y标签text(length(xx),xx(end),num2str(xx(end)),'fontsize',fs)%显示结果［图示］用下标作为自变量画迭代的折线。如P0_20_1图所示，当最大误差为10-4时，需要迭代19次才能达到精度，超越方程的解为27.539。［算法］方法二：用求零函数和求解函数。将方程改为函数f(x)=2ln(x)-3-100xMATLAB求零函数为fzero，fzero函数的格式之一是x=fzero(f,x0)其中，f表示求解的函数文件，x0是估计值。fzero函数的格式之二是x=fzero(f,［x1,x2］)

其中，x1和x2表示零点的范围。另外MATLAB还有求解函数solve，［程序］P1_2fzero.m如下。%超越方程的求法clearx=10:0.1:100;f=inline('2*log(x)-3T00./x')figureplot(x,f(x),'LineWidth',2)gridonx0=fzero(f,[20,30]);%x0=fzero(f,20);holdonplot(x0,f(x0),'.')计算非线性方程和方程组的符号解。%清除变量%自变量向量%定义内线函数用的是字符窜%创建图形窗口%画曲线%加网格%求方程的零点%求方程的零点%保持图像%画零点%标题%x标签%y标签%标记零点%求超越方程的符号解title('超越方程的解',计算非线性方程和方程组的符号解。%清除变量%自变量向量%定义内线函数用的是字符窜%创建图形窗口%画曲线%加网格%求方程的零点%求方程的零点%保持图像%画零点%标题%x标签%y标签%标记零点%求超越方程的符号解%清除变量%清除变量%间隔%自变量向量%原函数%通过差分求导数导数的计算正弦函数J=siM的导数是余弦函数V=cosx，余弦函数的导数是负的正弦函数，用MATLAB的数值导数和符号导数求正弦函数的一阶和二阶导数，并与其解析解进行比较。［程序］P2diff.m如下。%正弦函数导数的计算方法cleardx=0.01*2*pi;x=0:dx:2*pi;y=sin(x);f1=diff(y)/dx;f1=［f1(1),(f1(1:end-1)+f1(2:end))/2,f1(end)］;%求平均值figure %创建图形窗口

plot(x,cos(x),x,f1,'.')%plot(x,cos(x),x(1:end-1),f1,'.')%plot(x,cos(x),x(2:end),f1,'.')symssxy=sin(sx);dy_dx=diff(y);df1=subs(dy_dx,sx,x);holdonplot(x,df1,'ro')%画一阶导数和数值差分曲线%数值导数(点)偏左%数值导数(点)偏右%定义符号变量%建立符号函数%求符号导数%符号替换数值%保持图像%画符号导数曲线%加网格gridonlegend('解析导数'，’数值差分'，’符号导数'%画一阶导数和数值差分曲线%数值导数(点)偏左%数值导数(点)偏右%定义符号变量%建立符号函数%求符号导数%符号替换数值%保持图像%画符号导数曲线%加网格f2=diff(f1)/dx; %通过差分求导数f2=［f2(1),(f2(1:end-1)+f2(2:end))/2,f2(end)］;%求平均值d2y_dx2=diff(y,2); %求二阶符号导数df2=subs(d2y_dx2,sx,x); %符号替换数值figure %创建图形窗口plot(x,-sin(x),x,f2,'.',x,df2,'o') %画二阶导数和差分以及符号导数曲线gridon %加网格legend('解析导数'，'数值差分'，'符号导数',4)%图例title('正弦函数的二阶导数','FontSize',16)%加标题［图示］(1)如P2a图所示，正弦函数的一阶导数的数值解(点)与解析解(线)符合得很好。P2a图 P2b图(2)如P2a图 P2b图积分的计算求证：函数j=eax,vnbx的积分为S=——— eax(asinbx-bcosbx)+Ca2+b2其中a=-0.5,b=2。积分下限为0。上限为x，画出定积分的函数曲线。［证明］利用分部积分得

S=feaxsinbxdx=-fsinbxdeax=—{eaxsinbx-bfeaxcosbxdx}

aa={eaxsinbx-—\cosbxdeax}=-!-{eaxsinbx-—[eaxcosbx+bfeaxsinbxdx]}aa aa□ 1 b2S-—eax(asinbx一bcosbx)一一S

a2 a2由此可证不定积分。当x由此可证不定积分。当x=0时C=——

a2+b2因此，从0开始的积分为eax(asinbx-beax(asinbx-bcosbx+b)a2+b2利用复数积分的方法更简单。由于feax+ibxdx= eax+ibx+C'= eax+ibx+C'a+ib a2+b2其中仗表示复常数。根据欧拉公式eix=cosx+isinx，上式两边取虚部即可证明同一结果。上式两边取实部还可证明1eaxcosbxdx= eax(bsinbx+acosbx)+Ca2+b2［算法］设被积函数为J=f(x)，取间隔为Ax，取上限为x=nAx，则积分可用求和公式近似表示S=乎f(x)Axii=1积分既能用上式近似计算，也能用积分的解析式计算，还能用数值积分和符号积分计算。［程序］P3quad.m如下。%数值积分和符号积分方法%清除变量%指数的常数%正弦函数的常数%清除变量%指数的常数%正弦函数的常数%间隔%上限%自变量向量b=2;dx=0.1;xm=6;x=0:dx:xm;s1=(exp(a*x).*(-b*cos(b*x)+a*sin(b*x))+b)/(a“2+b“2);%积分的解析解y=exp(a*x).*sin(b*x); %被积函数s2=cumtrapz(y)*dx; %梯形法积分figure %创建图形窗口plot(x,s1,x,s2,'.') %画积分曲线gridon %加网格s=['exp(',num2str(a),'*x).*sin(',num2str(b),'*x)'];%被积分函数字符串f=inline(s); %化为内线函数，才可以被调用(画成)s3=0; %第1个积分值fori=2:length(x) %按自变量循环

s3=[s3,quad(f,0,x(i))];endholdonplot(x,s3,'or')s3=[s3,quad(f,0,x(i))];endholdonplot(x,s3,'or')symssasbsxss=exp(sa*sx)*sin(sb*sx);sy=int(ss,sx)ssy=subs(sy,{sa,sb},{a,b});s4=subs(ssy,sx,x);%结束循环%保持图像%画数值积分曲线%定义符号变量%被积符号函数%对sx进行符号积分%替换常数%替换向量因为sx与sa，sb的长度不一样，不能同时替代plot(x,s4-s4(1),'ko','MarkerSize',10)%画符号积分曲线tit=['\ity\rm=e"{',num2str(a),'}\it"x\rmsin',num2str(b),'\itx'];%形成数学公式：''表示上标title(［tit,'\rm的积分'］,'FontSize',16)%标题legend('公式法'，’梯形法'，’数值法'，’符号法',4)%加图例(4表示右下角，0表示电脑选择最佳位置，-1把图例放到外面)［图示］如P3图所示，梯形法积分(点)与积分的解析解(线)符合得很好，微分方程的求解方法(1)求一阶微分方程的解dy=2ydxx+1当x=0时，y=2，这是初始条件。用微分方程的数值解和符号解画出函数曲线，并与解析解进行比较。(2)求二阶微分方程的解d2y 2xdy- =0dx2 x2+1dxP3图初始条件为y(0)=1，y'(0)=2。用微分方程的数值解和符号解画出函数曲线和导数的曲线，并与解析解进行比较。P3图［解析］(1)分离变量得dy_2dx

yx+1积分得lny=2ln(x+1)+C利用初始条件可得C=ln2，因此y=2(x+1)2［程序］P4_1ode.m如下。%一阶常微分方程的解析解，数值解和符号解clear %清除变量x=linspace(0,2,50); %自变量向量y1=2*(x+1)."2; %解析解

f=inline('2*y/(x+1)'); %微分方程右边化为内线函数［x2,y2］=ode45(f,x,2); %求微分方程的数值解(ode常微分方程)ys=dsolve('Dy-2*y/(x+1)','y(0)=2','x')%求微分方程的符号特解dsolve(初始条件决定)Dy=dy_dxy3=subs(ys,'x',x);figureplot(x,y1,x,y2,'.',x,y3,'o')gridonlegend('解析解'，'数值解'，'符号解',4)xlabel('\itx','FontSize',16)%将符号改为向量求数值解%创建图形窗口%画曲线%加网格%图例%横坐标%纵坐标ylabel('\ity','FontSize',16)title(%将符号改为向量求数值解%创建图形窗口%画曲线%加网格%图例%横坐标%纵坐标［图示］如P4_1图所示，一阶微分方程的数值解(点)和符号解(圈)与解析解(线)符合得很好。P4_1图 P4_2图［解析］(2)由于y，=dy/dx，分离变量得dyf2xdx — y'x2+1积分得lny'-ln(x2+1)=C1当x=0时，y，=2，所以C1=ln2，因此y'=2(x2+1)再积分 2 -一y=1(2x2+2)dx=3x3+2x+C当x=0时，y=1，所以。2=1，因此y=—x3+2x+13设y(1)=y，y(2)=dy/dx，可得两个一阶微分方程dy⑴_6 dy⑵_2xy(2)云=y(2)，=日?将两个一阶微分方程设计成函数文件，以便求数值解。［程序］P0_23_2ode.m如下。

%二阶常微分方程的解析解，数值解和符号解clear %清除变量x=linspace(0,3,30); %自变量向量y1=1+2*x+2*x.”3/3; %解析解dy1=2*x.”2+2; %解析解的导数[x2,Y]=ode45('p4_2fun',x,[1,2]); %求微分方程的数值解[1,2]1，表示第一个初始条件默认了在(0，0)点附初值y2=Y(:,1); %取出函数‘：，1’表示第一列所有元素取出dy2=Y(:,2); %取出导数%符号导数%将符号改为向量求函数的数值解%将符号改为向量求导数的数值解%创建图形窗口%选子图%画曲线%加网格%图例%横坐标%纵坐标ys=dsolve('D2y-2*x*Dy/(x”2+1)','y(0)=1','Dy(0)=2','x')%求微分方程的符号特解((由初始条件决定)D2y标对y求二阶倒数’x'指出求得自变量是xdy=diff(ys);y3=subs(ys,'x',x);dy3=subs(dy,'x',x);figuresubplot(2,1,1)plot(x,y1,x2,y2,'.',x,y3,'o')gridonlegend('解析解'，'数值解'，'符号解',2)xlabel(%符号导数%将符号改为向量求函数的数值解%将符号改为向量求导数的数值解%创建图形窗口%选子图%画曲线%加网格%图例%横坐标%纵坐标title('二阶常微分方程解的函数','FontSize',16)%标题subplot(2,1,2) %选子图plot(x,dy1,x2,dy2,'.',x,dy3,'o') %画曲线gridon %加网格legend('解析解'，'数值解'，'符号解',2) %图例title('二阶常微分方程解的导数','FontSize',16)%标题xlabel('\itx','FontSize',16) %横坐标ylabel('d{\ity}/d\itx','FontSize',16) %纵坐标程序在求微分方程的数值解时将调用函数文件P4_2fun.m%二阶常微分方程的数值解的函数文件functionf=fun(x,y)f=[y(2); %一阶导数表达式2*x*y(2)/(x"2+1)]; %二阶导数表达式[图示]如P4_2图所示，二阶微分方程的数值解(包括函数和导数)和符号解与解析解都符合得很好。偏导数的计算和等量异号点电荷的电场两个异号点电荷带电量为±0。＞0)，相距为2s画出电场线和等势线。[解析]如B5图所示，等量异号点电荷在场点PXy)产生的电势为其中，k为静电力常量，。和r2是场点P到电荷的距离(2)电场强度可根据电势梯度计算其中，劈形算符为0=—i+—

dxE=-▽Uj+1kdy dz在xy平面上，场强只有两个分量E=-竺x dx竺dy两个点电荷在P点产生的电场强度的大小分别为场强的两个分量也能根据公式计算kQE=Ecos6-E场强的两个分量也能根据公式计算kQE=Ecos6-Ecos6=住")-kQ(x-a)

x1122r3 r3i2E=Esin6-Esin6=噂-您

y1 1 2 2r3 r312［算法］取a为坐标单位，做无量纲处理。则电势可表示为其中，U0=kQ/a。U0是Q在原点产生的电势r*=£=J(x*+1)2+y*2厂…1 1、U=Uo(r*-r*)12作为电势的单位。rr*=r「和r2*是约化距离(6)(7a)(7b)(1*)(2*)其中，x*=x/a，y*=y/a。x*和y*是无量纲的坐标或约化坐标。场强的x分量用梯度可表示为E=dU=-%dU*

x dx ad(x/a)(5a*)E=-E吐(5a*)x 0dx*其中，E0=U0/a，U*=U/Uo。Eo是场强的单位，U*是无量纲的电势。同理可得dU*(5b*)=-E (5b*)0dy*两个点电荷的电场强度的两个分量用公式可表示为X*+1X*-1 y*)*气-E°(~rr—"TT)，y一气(萩—元)12r*3ir*32(7*)%清除变量%横坐标范围%清除变量%横坐标范围%纵坐标范围%横坐标向量%纵坐标向量%坐标网点（矩阵）%左边第一个正电荷到场点的距离%右边第二个负电荷到场点的距离%计算电势%等势线的电势向量%创建图形窗口%画等势线并取等势线的坐标%标记等势线的值%加网格%保持图像%画水平和竖直线%[Ex,Ey]=gradient(-U);dth=20;th=(dth:dth:360-dth)*pi/180;r0=0.1;x0=r0*cos(th);y0=r0*sin(th);streamline(X,Y,Ex,Ey,x0-1,y0)streamline(X,Y,-Ex,-Ey,x0+1,y0)荷做正电荷处理)将物理量无量纲化之后，只要作纯数值计算就行了。MATLAB的梯度函数gradient可直接计算场强的数值分量，场强的数值解和解析解可相互比较。等势线可根据等值线指令contour绘制，电场线可根据流线指令streamline绘制。［程序］P0_24gradient.m如下。%等量异号点电荷的电场线和等势线（请在“创建图形窗口”处设置断点，以观察画图过程）clearxm=2.5;ym=2;x=linspace(-xm,xm,400);y=linspace(-ym,ym,400);[X,Y]=meshgrid(x,y);R1=sqrt((X+1)."2+Y."2);R2=sqrt((X-1)."2+Y."2);U=1./R1-1./R2