[程序]P1_1abs.m如下。
%超越方程的迭代算法
clear%清除变量
x0=30;%初始值
xx=[];%空向量
while1%无限循环
x=100/（2*log（x0）-3）;%迭代运算
xx=[xx,x];%连接结果
iflength（xx）>1000,break,end%如果项数太多则退出循环（暗示发散）
ifabs（x0-x）<1e-4,break,end%当精度足够高时退出循环
x0=x;%替换初值
end%结束循环
figure%创建图形窗口
plot（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。
[算法]方法二：
用求零函数和求解函数。
将方程改为函数
MATLAB求零函数为fzero，fzero函数的格式之一是
x=fzero（f,x0）
其中，f表示求解的函数文件，x0是估计值。
fzero函数的格式之二是
x=fzero（f,[x1,x2]）
其中，x1和x2表示零点的范围。
另外MATLAB还有求解函数solve，计算非线性方程和方程组的符号解。
[程序]P1_2fzero.m如下。
%超越方程的求法
clear%清除变量
x=10:
0.1:
100;%自变量向量
f=inline（'2*log（x）-3-100./x'）%定义内线函数用的是字符窜
figure%创建图形窗口
plot（x,f（x）,'LineWidth',2）%画曲线
gridon%加网格
x0=fzero（f,[20,30]）;%求方程的零点
%x0=fzero（f,20）;%求方程的零点
holdon%保持图像
plot（x0,f（x0）,'.'）%画零点
title（'超越方程的解','fontsize',16）%标题
xlabel（'\itx','fontsize',16）%x标签
ylabel（'\itf','fontsize',16）%y标签
text（x0,0,num2str（x0）,'fontsize',16）%标记零点
x0=solve（'2*log（x）-3-100./x'）%求超越方程的符号解
plot（double（x0）,0,'o'）%再画零点（double是把字符转化成数字）
P1_1图P1_2图
2.导数的计算
正弦函数y=sinx的导数是余弦函数y'=cosx，余弦函数的导数是负的正弦函数，用MATLAB的数值导数和符号导数求正弦函数的一阶和二阶导数，并与其解析解进行比较。
[程序]P2diff.m如下。
%正弦函数导数的计算方法
clear%清除变量
dx=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,'.'）%数值导数（点）偏右
symssx%定义符号变量
y=sin（sx）;%建立符号函数
dy_dx=diff（y）;%求符号导数
df1=subs（dy_dx,sx,x）;%符号替换数值
holdon%保持图像
plot（x,df1,'ro'）%画符号导数曲线
gridon%加网格
legend（'解析导数','数值差分','符号导数',4）%图例
title（'正弦函数的一阶导数','FontSize',16）%加标题
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图所示，正弦函数的一阶导数的数值解（点）与解析解（线）符合得很好。
（2）如P2b图所示，正弦函数的二阶导数的数值解（点）和符号解（圈）与解析解（线）符合得很好，不过二阶数值导数在端点与精确值有一点偏离。
P2a图P2b图
3.积分的计算
求证：
函数y=eaxsinbx的积分为
其中a=-0.5，b=2。
积分下限为0。
上限为x，画出定积分的函数曲线。
[证明]利用分部积分得
即
由此可证不定积分。
当x=0时，S应该为零，所以
因此，从0开始的积分为
利用复数积分的方法更简单。
由于
其中C'表示复常数。
根据欧拉公式eix=cosx+isinx，上式两边取虚部即可证明同一结果。
上式两边取实部还可证明
[算法]设被积函数为y=f（x），取间隔为Δx，取上限为x=nΔx，则积分可用求和公式近似表示
积分既能用上式近似计算，也能用积分的解析式计算，还能用数值积分和符号积分计算。
[程序]P3quad.m如下。
%数值积分和符号积分方法
clear%清除变量
a=-0.5;%指数的常数
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））];%连接积分quad对f积分下限0，上限x（i）
end%结束循环
holdon%保持图像
plot（x,s3,'or'）%画数值积分曲线
symssasbsx%定义符号变量
ss=exp（sa*sx）*sin（sb*sx）;%被积符号函数
sy=int（ss,sx）%对sx进行符号积分
ssy=subs（sy,{sa,sb},{a,b}）;%替换常数
s4=subs（ssy,sx,x）;%替换向量因为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图所示，梯形法积分（点）与积分的解析解（线）符合得很好，
4.微分方程的求解方法
（1）求一阶微分方程的解
当x=0时，y=2，这是初始条件。
用微分方程的数值解和符号解画出函数曲线，并与解析解进行比较。
（2）求二阶微分方程的解
P3图
初始条件为y（0）=1，y'（0）=2。
用微分方程的数值解和符号解画出函数曲线和导数的曲线，并与解析解进行比较。
[解析]
（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_dx
y3=subs（ys,'x',x）;%将符号改为向量求数值解
figure%创建图形窗口
plot（x,y1,x,y2,'.',x,y3,'o'）%画曲线
gridon%加网格
legend（'解析解','数值解','符号解',4）%图例
xlabel（'\itx','FontSize',16）%横坐标
ylabel（'\ity','FontSize',16）%纵坐标
title（'一阶常微分方程的解','FontSize',16）%标题
[图示]如P4_1图所示，一阶微分方程的数值解（点）和符号解（圈）与解析解（线）符合得很好。
P4_1图P4_2图
[解析]
（2）由于y'=dy/dx，分离变量得
积分得
lny'-ln（x2+1）=C1
当x=0时，y'=2，所以C1=ln2，因此
y'=2（x2+1）
再积分
当x=0时，y=1，所以C2=1，因此
设y
（1）=y，y
（2）=dy/dx，可得两个一阶微分方程
，
将两个一阶微分方程设计成函数文件，以便求数值解。
[程序]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'指出求得自变量是x
dy=diff（ys）;%符号导数
y3=subs（ys,'x',x）;%将符号改为向量求函数的数值解
dy3=subs（dy,'x',x）;%将符号改为向量求导数的数值解
figure%创建图形窗口
subplot（2,1,1）%选子图
plot（x,y1,x2,y2,'.',x,y3,'o'）%画曲线
gridon%加网格
legend（'解析解','数值解','符号解',2）%图例
xlabel（'\itx','FontSize',16）%横坐标
ylabel（'\ity','FontSize',16）%纵坐标
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图所示，二阶微分方程的数值解（包括函数和导数）和符号解与解析解都符合得很好。
5.偏导数的计算和等量异号点电荷的电场
两个异号点电荷带电量为±Q（Q>0），相距为2a，画出电场线和等势线。
[解析]如B5图所示，等量异号点电荷在场点P（x,y）产生的电势为
（1）
其中，k为静电力常量，r1和r2是场点P到电荷的距离
，
（2）
电场强度可根据电势梯度计算
E=-▽U（3）
其中，劈形算符为
（4）
在xy平面上，场强只有两个分量
，
（5）
两个点电荷在P点产生的电场强度的大小分别为
，
（6）
场强的两个分量也能根据公式计算
（7a）
（7b）
[算法]取a为坐标单位，做无量纲处理。
则电势可表示为
（1*）
其中，U0=kQ/a。
U0是Q在原点产生的电势，作为电势的单位。
r1*和r2*是约化距离
，
（2*）
其中，x*=x/a，y*=y/a。
x*和y*是无量纲的坐标或约化坐标。
场强的x分量用梯度可表示为
即
（5a*）
其中，E0=U0/a，U*=U/U0。
E0是场强的单位，U*是无量纲的电势。
同理可得
（5b*）
两个点电荷的电场强度的两个分量用公式可表示为
，
（7*）
将物理量无量纲化之后，只要作纯数值计算就行了。
MATLAB的梯度函数gradient可直接计算场强的数值分量，场强的数值解和解析解可相互比较。
等势线可根据等值线指令contour绘制，电场线可根据流线指令streamline绘制。
[程序]P0_24gradient.m如下。
%等量异号点电荷的电场线和等势线（请在“创建图形窗口”处设置断点，以观察画图过程）
clear%清除变量
xm=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;%计算电势
u=-4:
0.5:
4;%等势线的电势向量
figure%创建图形窗口
C=contour（X,Y,U,u,'LineWidth',2）;%画等势线并取等势线的坐标
clabel（C,'fontsize',16）%标记等势线的值
gridon%加网格
holdon%保持图像
plot（[-xm;xm],[0;0],[0;0],[-ym;ym]）%画水平和竖直线
（一）安全评价的内涵plot（-1,0,'o',1,0,'o','MarkerSize',12）%画电荷
1）按类型分。
环境标准按类型分为环境质量标准、污染物排放标准（或控制标准）、环境基础标准、环境检测方法标准、环境标准样品标准。
[Ex,Ey]=gradient（-U,x
（2）-x
（1）,y
（2）-y
（1））;%用电势梯度求场强的两个分量
2.早期介入原则；%[Ex,Ey]=gradient（-U）;%用电势梯度求场强的两个分量
（5）阐述划分评价单元的原则、分析过程等。
dth=20;%电场线角度间隔
（1）安全预评价。
th=（dth:
dth:
360-dth）*pi/180;%电场线的起始角度
r0=0.1;%电场线起点半径
x0=r0*cos（th）;%电场线的起点横坐标
y0=r0*sin（th）;%电场线的起点纵坐标
1.规划环境影响评价的报审streamline（X,Y,Ex,Ey,x0-1,y0）%画左边电场线（中间部分达到右边）
streamline（X,Y,-Ex,-Ey,x0+1,y0）%画右边电场线（中间部分达到左边）（-Ex是把负电荷做正电荷处理）
4.建设项目环境影响评价文件的分级审批axisequaltight%使坐标间隔相等
（4）预防或者减轻不良环境影响的对策和措施的合理性和有效性；title（'等量异号点电荷的电场线和等势线','fontsize',20）%显示标题
xlabel（'\itx/a','fontsize',16）%显示横坐标
ylabel（'\ity/a','fontsize',16）%显示纵坐标
text（-xm,ym-0.5,'电势单位:
\itkQ/a','fontsize',16）%显示电势单位
安全评价是落实“安全第一，预防为主，综合治理”方针的重要技术保障，是安全生产监督管理的重要手段。
Ex=（X+1）./R1.^3-（X-1）./R2.^3;%用公式求场强的x分量
Ey=Y./R1.^3-Y./R2.^3;%用公式求场强的y分量
streamline（X,Y,Ex,Ey,x0-1,y0）%重画左边电场线（曲线重合）
streamline（X,Y,-Ex,-Ey,x0+1,y0）%重画右上电场线（曲线重合）
1.依法评价原则；[图示]如P0_24图所示，左边表示正电荷，右边表示负电荷，等量异号点电荷的电场线和等势线关于原点是对称分布的。
电场线从正电荷出发，终止于负电荷。
电场线与等势线垂直，任何两条电场线都不相交。
除了电势为零的直线外，等势线分别包围着各自的电荷。
电场强度大的地方，电场线较密，等势线也较密。
升级会员 