MATLAB数学建模6常用计算方法文档格式.docx
《MATLAB数学建模6常用计算方法文档格式.docx》由会员分享，可在线阅读，更多相关《MATLAB数学建模6常用计算方法文档格式.docx（15页珍藏版）》请在冰豆网上搜索。
title（'
超越方程的迭代折线'
FontSize'
fs）%标题
xlabel（'
\itn'
fs）%x标签
ylabel（'
\itx'
fs）%y标签
text（length（xx）,xx（end）,num2str（xx（end））,'
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如下。
%超越方程的求法
x=10:
0.1:
100;
%自变量向量
f=inline（'
2*log（x）-3-100./x'
）%定义内线函数
plot（x,f（x）,'
2）%画曲线
x0=fzero（f,[20,30]）;
%求方程的零点
%x0=fzero（f,20）;
holdon%保持图像
plot（x0,f（x0）,'
.'
）%画零点
超越方程的解'
16）%标题
16）%x标签
\itf'
16）%y标签
text（x0,0,num2str（x0）,'
16）%标记零点
x0=solve（'
）%求超越方程的符号解
plot（double（x0）,0,'
o'
）%再画零点
P1_1图P1_2图
2.导数的计算
正弦函数y=sinx的导数是余弦函数y'
=cosx，余弦函数的导数是负的正弦函数，用MATLAB的数值导数和符号导数求正弦函数的一阶和二阶导数，并与其解析解进行比较。
[程序]P2diff.m如下。
%正弦函数导数的计算方法
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）];
%求平均值
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）;
%符号替换数值
plot（x,df1,'
ro'
）%画符号导数曲线
legend（'
解析导数'
数值差分'
符号导数'
4）%图例
正弦函数的一阶导数'
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）;
plot（x,-sin（x）,x,f2,'
x,df2,'
）%画二阶导数和差分以及符号导数曲线
正弦函数的二阶导数'
[图示]
（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如下。
%数值积分和符号积分方法
a=-0.5;
%指数的常数
b=2;
%正弦函数的常数
dx=0.1;
xm=6;
%上限
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;
%梯形法积分
plot（x,s1,x,s2,'
）%画积分曲线
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））];
%连接积分
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）;
%替换向量
plot（x,s4-s4
（1）,'
ko'
10）%画符号积分曲线
tit=['
\ity\rm=e^{'
}\it^x\rmsin'
%公式
title（[tit,'
\rm的积分'
],'
公式法'
梯形法'
数值法'
符号法'
4）%加图例
[图示]如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如下。
%一阶常微分方程的解析解,数值解和符号解
x=linspace（0,2,50）;
y1=2*（x+1）.^2;
%解析解
2*y/（x+1）'
）;
%微分方程右边化为内线函数
[x2,y2]=ode45（f,x,2）;
%求微分方程的数值解
ys=dsolve（'
Dy-2*y/（x+1）'
y（0）=2'
x'
）%求微分方程的符号特解（由初始条件决定）
y3=subs（ys,'
x）;
%将符号改为向量求数值解
plot（x,y1,x,y2,'
x,y3,'
）%画曲线
解析解'
数值解'
符号解'
4）%图例
16）%横坐标
\ity'
16）%纵坐标
一阶常微分方程的解'
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如下。
%二阶常微分方程的解析解,数值解和符号解
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]）;
y2=Y（:
1）;
%取出函数
dy2=Y（:
2）;
%取出导数
D2y-2*x*Dy/（x^2+1）'
y（0）=1'
Dy（0）=2'
）%求微分方程的符号特解（（由初始条件决定）
dy=diff（ys）;
%符号导数
%将符号改为向量求函数的数值解
dy3=subs（dy,'
%将符号改为向量求导数的数值解
subplot（2,1,1）%选子图
plot（x,y1,x2,y2,'
2）%图例
二阶常微分方程解的函数'
subplot（2,1,2）%选子图
plot（x,dy1,x2,dy2,'
x,dy3,'
二阶常微分方程解的导数'
d{\ity}/d\itx'
程序在求微分方程的数值解时将调用函数文件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如下。
%等量异号点电荷的电场线和等势线（请在“创建图形窗口”处设置断点，以观察画图过程）
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;
%等势线的电势向量
C=contour（X,Y,U,u,'
%画等势线并取等势线的坐标
clabel（C,'
16）%标记等势线的值
plot（[-xm;
xm],[0;
0],[0;
0],[-ym;
ym]）%画水平和竖直线
plot（-1,0,'
1,0,'
12）%画电荷
[Ex,Ey]=gradient（-U,x
（2）-x
（1）,y
（2）-y
（1））;
%用电势梯度求场强的两个分量
%[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）%画右边电场线（中间部分达到左边）
axisequaltight%使坐标间隔相等
等量异号点电荷的电场线和等势线'
20）%显示标题
\itx/a'
16）%显示横坐标
\ity/a'
16）%显示纵坐标
text（-xm,ym-0.5,'
电势单位:
\itkQ/a'
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）%重画右上电场线（曲线重合）
[图示]如P0_24图所示，左边表示正电荷，右边表示负电荷，等量异号点电荷的电场线和等势线关于原点是对称分布的。
电场线从正电荷出发，终止于负电荷。
电场线与等势线垂直，任何两条电场线都不相交。
除了电势为零的直线外，等势线分别包围着各自的电荷。
电场强度大的地方，电场线较密，等势线也较密。
升级会员 