常用計算方法1.超越方程的求解一超越方程為x(2lnx-3)-100=0求超越方程的解。［算法］方法一：用迭代算法。將方程改為100X= 2ln(X0)-3其中x0是一個初始值，由此計算終值x。取最大誤差為e=10-4，當(dāng)|x-x0l>e時，就用x的值換成x0的值，重新進行計算；否則lx-x0l<e為止。［程序］P1_1abs.m如下。%超越方程的迭代算法%清除變量%初始值%空向量%清除變量%初始值%空向量%無限循環(huán)%迭代運算%連接結(jié)果%如果項數(shù)太多則退出循環(huán)(暗示發(fā)散)%當(dāng)精度足夠高時退出循環(huán)%替換初值%結(jié)束循環(huán)%創(chuàng)建圖形窗口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 %加網(wǎng)格fs=16; %字體大小title('超越方程的迭代折線’,'fontsize',fs)%標(biāo)題xlabel('\itn', 'fontsize',fs) %x標(biāo)簽ylabel('\itx', 'fontsize',fs) %y標(biāo)簽text(length(xx),xx(end),num2str(xx(end)),'fontsize',fs)%顯示結(jié)果［圖示］用下標(biāo)作為自變量畫迭代的折線。如P0_20_1圖所示，當(dāng)最大誤差為10-4時，需要迭代19次才能達到精度，超越方程的解為27.539。［算法］方法二：用求零函數(shù)和求解函數(shù)。將方程改為函數(shù)f(x)=2ln(x)-3-100xMATLAB求零函數(shù)為fzero，fzero函數(shù)的格式之一是x=fzero(f,x0)其中，f表示求解的函數(shù)文件，x0是估計值。fzero函數(shù)的格式之二是x=fzero(f,［x1,x2］)

其中，x1和x2表示零點的范圍。另外MATLAB還有求解函數(shù)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),'.')計算非線性方程和方程組的符號解。%清除變量%自變量向量%定義內(nèi)線函數(shù)用的是字符竄%創(chuàng)建圖形窗口%畫曲線%加網(wǎng)格%求方程的零點%求方程的零點%保持圖像%畫零點%標(biāo)題%x標(biāo)簽%y標(biāo)簽%標(biāo)記零點%求超越方程的符號解title('超越方程的解',計算非線性方程和方程組的符號解。%清除變量%自變量向量%定義內(nèi)線函數(shù)用的是字符竄%創(chuàng)建圖形窗口%畫曲線%加網(wǎng)格%求方程的零點%求方程的零點%保持圖像%畫零點%標(biāo)題%x標(biāo)簽%y標(biāo)簽%標(biāo)記零點%求超越方程的符號解%清除變量%清除變量%間隔%自變量向量%原函數(shù)%通過差分求導(dǎo)數(shù)導(dǎo)數(shù)的計算正弦函數(shù)J=siM的導(dǎo)數(shù)是余弦函數(shù)V=cosx，余弦函數(shù)的導(dǎo)數(shù)是負的正弦函數(shù)，用MATLAB的數(shù)值導(dǎo)數(shù)和符號導(dǎo)數(shù)求正弦函數(shù)的一階和二階導(dǎo)數(shù)，并與其解析解進行比較。［程序］P2diff.m如下。%正弦函數(shù)導(dǎo)數(shù)的計算方法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 %創(chuàng)建圖形窗口

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')%畫一階導(dǎo)數(shù)和數(shù)值差分曲線%數(shù)值導(dǎo)數(shù)(點)偏左%數(shù)值導(dǎo)數(shù)(點)偏右%定義符號變量%建立符號函數(shù)%求符號導(dǎo)數(shù)%符號替換數(shù)值%保持圖像%畫符號導(dǎo)數(shù)曲線%加網(wǎng)格gridonlegend('解析導(dǎo)數(shù)'，’數(shù)值差分'，’符號導(dǎo)數(shù)'%畫一階導(dǎo)數(shù)和數(shù)值差分曲線%數(shù)值導(dǎo)數(shù)(點)偏左%數(shù)值導(dǎo)數(shù)(點)偏右%定義符號變量%建立符號函數(shù)%求符號導(dǎo)數(shù)%符號替換數(shù)值%保持圖像%畫符號導(dǎo)數(shù)曲線%加網(wǎng)格f2=diff(f1)/dx; %通過差分求導(dǎo)數(shù)f2=［f2(1),(f2(1:end-1)+f2(2:end))/2,f2(end)］;%求平均值d2y_dx2=diff(y,2); %求二階符號導(dǎo)數(shù)df2=subs(d2y_dx2,sx,x); %符號替換數(shù)值figure %創(chuàng)建圖形窗口plot(x,-sin(x),x,f2,'.',x,df2,'o') %畫二階導(dǎo)數(shù)和差分以及符號導(dǎo)數(shù)曲線gridon %加網(wǎng)格legend('解析導(dǎo)數(shù)'，'數(shù)值差分'，'符號導(dǎo)數(shù)',4)%圖例title('正弦函數(shù)的二階導(dǎo)數(shù)','FontSize',16)%加標(biāo)題［圖示］(1)如P2a圖所示，正弦函數(shù)的一階導(dǎo)數(shù)的數(shù)值解(點)與解析解(線)符合得很好。P2a圖 P2b圖(2)如P2a圖 P2b圖積分的計算求證：函數(shù)j=eax,vnbx的積分為S=——— eax(asinbx-bcosbx)+Ca2+b2其中a=-0.5,b=2。積分下限為0。上限為x，畫出定積分的函數(shù)曲線。［證明］利用分部積分得

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

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

a2 a2由此可證不定積分。當(dāng)x由此可證不定積分。當(dāng)x=0時C=——

a2+b2因此，從0開始的積分為eax(asinbx-beax(asinbx-bcosbx+b)a2+b2利用復(fù)數(shù)積分的方法更簡單。由于feax+ibxdx= eax+ibx+C'= eax+ibx+C'a+ib a2+b2其中仗表示復(fù)常數(shù)。根據(jù)歐拉公式eix=cosx+isinx，上式兩邊取虛部即可證明同一結(jié)果。上式兩邊取實部還可證明1eaxcosbxdx= eax(bsinbx+acosbx)+Ca2+b2［算法］設(shè)被積函數(shù)為J=f(x)，取間隔為Ax，取上限為x=nAx，則積分可用求和公式近似表示S=乎f(x)Axii=1積分既能用上式近似計算，也能用積分的解析式計算，還能用數(shù)值積分和符號積分計算。［程序］P3quad.m如下。%數(shù)值積分和符號積分方法%清除變量%指數(shù)的常數(shù)%正弦函數(shù)的常數(shù)%清除變量%指數(shù)的常數(shù)%正弦函數(shù)的常數(shù)%間隔%上限%自變量向量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); %被積函數(shù)s2=cumtrapz(y)*dx; %梯形法積分figure %創(chuàng)建圖形窗口plot(x,s1,x,s2,'.') %畫積分曲線gridon %加網(wǎng)格s=['exp(',num2str(a),'*x).*sin(',num2str(b),'*x)'];%被積分函數(shù)字符串f=inline(s); %化為內(nèi)線函數(shù)，才可以被調(diào)用(畫成)s3=0; %第1個積分值fori=2:length(x) %按自變量循環(huán)

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);%結(jié)束循環(huán)%保持圖像%畫數(shù)值積分曲線%定義符號變量%被積符號函數(shù)%對sx進行符號積分%替換常數(shù)%替換向量因為sx與sa，sb的長度不一樣，不能同時替代plot(x,s4-s4(1),'ko','MarkerSize',10)%畫符號積分曲線tit=['\ity\rm=e"{',num2str(a),'}\it"x\rmsin',num2str(b),'\itx'];%形成數(shù)學(xué)公式：''表示上標(biāo)title(［tit,'\rm的積分'］,'FontSize',16)%標(biāo)題legend('公式法'，’梯形法'，’數(shù)值法'，’符號法',4)%加圖例(4表示右下角，0表示電腦選擇最佳位置，-1把圖例放到外面)［圖示］如P3圖所示，梯形法積分(點)與積分的解析解(線)符合得很好，微分方程的求解方法(1)求一階微分方程的解dy=2ydxx+1當(dāng)x=0時，y=2，這是初始條件。用微分方程的數(shù)值解和符號解畫出函數(shù)曲線，并與解析解進行比較。(2)求二階微分方程的解d2y 2xdy- =0dx2 x2+1dxP3圖初始條件為y(0)=1，y'(0)=2。用微分方程的數(shù)值解和符號解畫出函數(shù)曲線和導(dǎo)數(shù)的曲線，并與解析解進行比較。P3圖［解析］(1)分離變量得dy_2dx

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

f=inline('2*y/(x+1)'); %微分方程右邊化為內(nèi)線函數(shù)［x2,y2］=ode45(f,x,2); %求微分方程的數(shù)值解(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('解析解'，'數(shù)值解'，'符號解',4)xlabel('\itx','FontSize',16)%將符號改為向量求數(shù)值解%創(chuàng)建圖形窗口%畫曲線%加網(wǎng)格%圖例%橫坐標(biāo)%縱坐標(biāo)ylabel('\ity','FontSize',16)title(%將符號改為向量求數(shù)值解%創(chuàng)建圖形窗口%畫曲線%加網(wǎng)格%圖例%橫坐標(biāo)%縱坐標(biāo)［圖示］如P4_1圖所示，一階微分方程的數(shù)值解(點)和符號解(圈)與解析解(線)符合得很好。P4_1圖 P4_2圖［解析］(2)由于y，=dy/dx，分離變量得dyf2xdx — y'x2+1積分得lny'-ln(x2+1)=C1當(dāng)x=0時，y，=2，所以C1=ln2，因此y'=2(x2+1)再積分 2 -一y=1(2x2+2)dx=3x3+2x+C當(dāng)x=0時，y=1，所以。2=1，因此y=—x3+2x+13設(shè)y(1)=y，y(2)=dy/dx，可得兩個一階微分方程dy⑴_6 dy⑵_2xy(2)云=y(2)，=日?將兩個一階微分方程設(shè)計成函數(shù)文件，以便求數(shù)值解。［程序］P0_23_2ode.m如下。

%二階常微分方程的解析解，數(shù)值解和符號解clear %清除變量x=linspace(0,3,30); %自變量向量y1=1+2*x+2*x.”3/3; %解析解dy1=2*x.”2+2; %解析解的導(dǎo)數(shù)[x2,Y]=ode45('p4_2fun',x,[1,2]); %求微分方程的數(shù)值解[1,2]1，表示第一個初始條件默認了在(0，0)點附初值y2=Y(:,1); %取出函數(shù)‘：，1’表示第一列所有元素取出dy2=Y(:,2); %取出導(dǎo)數(shù)%符號導(dǎo)數(shù)%將符號改為向量求函數(shù)的數(shù)值解%將符號改為向量求導(dǎo)數(shù)的數(shù)值解%創(chuàng)建圖形窗口%選子圖%畫曲線%加網(wǎng)格%圖例%橫坐標(biāo)%縱坐標(biāo)ys=dsolve('D2y-2*x*Dy/(x”2+1)','y(0)=1','Dy(0)=2','x')%求微分方程的符號特解((由初始條件決定)D2y標(biāo)對y求二階倒數(shù)’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('解析解'，'數(shù)值解'，'符號解',2)xlabel(%符號導(dǎo)數(shù)%將符號改為向量求函數(shù)的數(shù)值解%將符號改為向量求導(dǎo)數(shù)的數(shù)值解%創(chuàng)建圖形窗口%選子圖%畫曲線%加網(wǎng)格%圖例%橫坐標(biāo)%縱坐標(biāo)title('二階常微分方程解的函數(shù)','FontSize',16)%標(biāo)題subplot(2,1,2) %選子圖plot(x,dy1,x2,dy2,'.',x,dy3,'o') %畫曲線gridon %加網(wǎng)格legend('解析解'，'數(shù)值解'，'符號解',2) %圖例title('二階常微分方程解的導(dǎo)數(shù)','FontSize',16)%標(biāo)題xlabel('\itx','FontSize',16) %橫坐標(biāo)ylabel('d{\ity}/d\itx','FontSize',16) %縱坐標(biāo)程序在求微分方程的數(shù)值解時將調(diào)用函數(shù)文件P4_2fun.m%二階常微分方程的數(shù)值解的函數(shù)文件functionf=fun(x,y)f=[y(2); %一階導(dǎo)數(shù)表達式2*x*y(2)/(x"2+1)]; %二階導(dǎo)數(shù)表達式[圖示]如P4_2圖所示，二階微分方程的數(shù)值解(包括函數(shù)和導(dǎo)數(shù))和符號解與解析解都符合得很好。偏導(dǎo)數(shù)的計算和等量異號點電荷的電場兩個異號點電荷帶電量為±0。＞0)，相距為2s畫出電場線和等勢線。[解析]如B5圖所示，等量異號點電荷在場點PXy)產(chǎn)生的電勢為其中，k為靜電力常量，。和r2是場點P到電荷的距離(2)電場強度可根據(jù)電勢梯度計算其中，劈形算符為0=—i+—

dxE=-▽Uj+1kdy dz在xy平面上，場強只有兩個分量E=-竺x dx竺dy兩個點電荷在P點產(chǎn)生的電場強度的大小分別為場強的兩個分量也能根據(jù)公式計算kQE=Ecos6-E場強的兩個分量也能根據(jù)公式計算kQE=Ecos6-Ecos6=住")-kQ(x-a)

x1122r3 r3i2E=Esin6-Esin6=噂-您

y1 1 2 2r3 r312［算法］取a為坐標(biāo)單位，做無量綱處理。則電勢可表示為其中，U0=kQ/a。U0是Q在原點產(chǎn)生的電勢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*是無量綱的坐標(biāo)或約化坐標(biāo)。場強的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*)%清除變量%橫坐標(biāo)范圍%清除變量%橫坐標(biāo)范圍%縱坐標(biāo)范圍%橫坐標(biāo)向量%縱坐標(biāo)向量%坐標(biāo)網(wǎng)點（矩陣）%左邊第一個正電荷到場點的距離%右邊第二個負電荷到場點的距離%計算電勢%等勢線的電勢向量%創(chuàng)建圖形窗口%畫等勢線并取等勢線的坐標(biāo)%標(biāo)記等勢線的值%加網(wǎng)格%保持圖像%畫水平和豎直線%[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)荷做正電荷處理)將物理量無量綱化之后，只要作純數(shù)值計算就行了。MATLAB的梯度函數(shù)gradient可直接計算場強的數(shù)值分量，場強的數(shù)值解和解析解可相互比較。等勢線可根據(jù)等值線指令contour繪制，電場線可根據(jù)流線指令streamline繪制。［程序］P0_24gradient.m如下。%等量異號點電荷的電場線和等勢線（請在“創(chuàng)建圖形窗口”處設(shè)置斷點，以觀察畫圖過程）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