美赛整理之Matlab的工程数学计算学习笔记（高等数学）

1.极限的定义和判别：
2.绘制特殊曲面
3.求两个空间曲面的交线
4.定积分的计算
5.多重积分的计算
- 1.截面法：
- 2.定义法
- - （1）先画图：
  - （2）分析：
- 3.蒙特卡洛法

1.极限的定义和判别：

由极限的定义可知：对于给定的f(x)f(x)f(x),∀ϵ>0,∃δ>0，\forall \epsilon>0,\exists \delta>0，∀ϵ>0,∃δ>0，使得：当∣x−x0∣<δ|x-x_0|<\delta∣x−x0∣<δ时，有∣f(x)−A∣<ϵ|f(x)-A|<\epsilon∣f(x)−A∣<ϵ成立，那么我们称：limx→x0f(x)=Alim_{x\rightarrow x_0}f(x)=Alimx→x0f(x)=A。

那么我们的关键就是要给出∀ϵ>0\forall \epsilon>0∀ϵ>0,我们都能求出对应的δ\deltaδ。

matlab的实现

clc,clear;
disp('判断A是否时f(x)的极限值\n');
A= input('请输入A的极限值:\n');
x_c = input('请输入对应的自变量的值：\n');
fxc = input('f(x)的表达式为，例如sin(x)/x\n','s');
delta = 1;n = 1;x = x_c - delta;
flag = 1;p = 1;
while flag == 1epsilon = input('请输入一个小鼠:\n');while abs(A - eval(fxc))>epsilon delta = delta/2;if p == 1 x = x_c - delta;elsex = x_c + delta;endif abs(delta)<epsdisp('没有找到这样的delta？');n = 0;breakendendif n == 0if p == 1disp('左极限不正确。');elsedisp('右极限不正确');endbreakendif n==1if p == 1disp('左极限可能正确');elsedisp('右极限可能正确');enddisp('delta的值是：');disp(delta);endflag = input('还要重新试一下吗？1：重新试一下/0：不试了\n');if flag == 1p = input('请告诉我这次要判断左极限还是右极限？1：左极限，0：右极限\n');end
end

2.绘制特殊曲面

若二次曲面的方程如下：
x2a2+y2b2+z2c2=d\frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2} = d a2x2+b2y2+c2z2=d
根据a,b,ca,b,ca,b,c的不同绘制出不同的图形。

在给定x,yx,yx,y要求求出zzz时，如果有开方的运算，会出现虚数，而且对于实数也会有两个符号不同的解。为了使得虚数不会出现在绘图中，我们把虚数换成·非数(NaNNaNNaN)

我们规定：如果一个z值只有只有这一点时虚数，我们就让它为0。如果这个数四周都是虚数，我们就让它时NaNNaNNaN

clc,clear;
a = input('a = ');
b = input('b = ');
c = input('c = ');
d = input('d = ');
N = input('划分网格点的数目为：');
xgrid = linspace(-abs(a),abs(a),N);
ygrid = linspace(-abs(b),abs(b),N);
[x,y] = meshgrid(xgrid,ygrid);
z = c*sqrt(d - y.^2/b^2-x.^2/a^2);
u=input('是否要画出全部图像？1：画出\0:不画\n');
z1 = real(z);
for i = 2:N-1for j = 2:N-1if(imag(z(i,j))~=0)z1(i,j) = 0;endif all(imag(z([i-1:i+1],[j-1:j+1]))~=0)z1(i,j) = NaN;endend
end
surf(x,y,z1);
hold on
if u == 1z2 = - z1;surf(x,y,z2);
end
axis([-abs(a),abs(a),-abs(b),abs(b),-abs(c),abs(c)]);
xlabel('x');ylabel('y');zlabel('z');
hold off

3.求两个空间曲面的交线

求以下两个方程的交线：
z=x2−2y2z=2x−3yz = x^2-2y^2\\ z = 2x-3y z=x2−2y2z=2x−3y

clc,clear;
[x,y] = meshgrid(-2:0.1:2);
z1 = x.^2 - 2*y.^2;
z2 = 2*x - 3*y;
mesh(x,y,z1);hold on;mesh(x,y,z2);
r0 = (abs(z1-z2)<0.01);
zz = r0.*z1;
yy = r0.*y;
xx = r0.*x;
plot3(xx(r0~=0),yy(r0~=0),zz(r0~=0),'--');

4.定积分的计算

先建立example.mexample.mexample.m函数。再调用Matlab的函数：s=quad(′example′,a,b)s = quad('example',a,b)s=quad(′example′,a,b)

5.多重积分的计算

1.截面法：

∫∫D(x2+y2)dσD={由x=1,y=x,y=0围成的区域}\int\int_D(x^2+y^2)d\sigma\\ D = \{由x = 1,y=x,y = 0围成的区域\} ∫∫D(x2+y2)dσD={由x=1,y=x,y=0围成的区域}

clc,clear;
fill([0,1,1,0],[0,0,1,0],'w');
hold on
fill([0.55,0.6,0.6,0.55,0.55],[0,0,0.6,0.55,0],'r');
dx = input('请告诉我输入步长dx = ');
dy = dx;
x = 0:dx:1;
lx = length(x);
f_value = @(x)(x);
for i = 1:lxx1 = (i-1)*dx;y1 = 0:dy:f_value(x1);ly = length(y1);f = (x1*ones(1,ly)).^2+y1.^2;s1(i) = trapz(f)*dy;
end
s = trapz(s1)*dx;
disp(s);

答案是1/31/31/3。

2.定义法

计算：
∫∫∫Ωxy2z3dxdydz\int\int\int_{\Omega}xy^2z^3dxdydz ∫∫∫Ωxy2z3dxdydz
其中：Ω\OmegaΩ是由x=1,y=x,z=xy,z=0x=1,y=x,z=xy,z=0x=1,y=x,z=xy,z=0围成的区域。

分析·：

（1）先画图：

[x,y] = meshgrid(0:0.05:1);
z1 = x.*y;                       %z1的表面
z2 = zeros(size(z1));            %z2的表面
mesh(x,y,z1);hold on
mesh(x,y,z2);
x1 = [0:0.02:1];
y1 = x1;
sx1 = length(x1);
zd = [zeros(1,sx1);x1.*y1];
line([x1;x1],[y1;y1],zd);
plot3([x1;x1],[y1;y1],zd,'*');
plot3(ones(2,sx1),[y1;y1],[zeros(1,sx1);y1],'o');
xlabel('x');
ylabel('y');
zlabel('z');

（2）分析：

∫∫∫Ωf(x,y,z)dV=limmax{ΔV}→0∑i=0n∑j=0m∑k=0tf(xi,yj,zk)dxidyjdzk\int\int\int_{\Omega}f(x,y,z)dV = lim_{max_{\{\Delta V\}}\rightarrow0}\sum_{i=0}^{n}\sum_{j=0}^{m}\sum_{k=0}^{t} f(x_i,y_j,z_k)dx_idy_jdz_k ∫∫∫Ωf(x,y,z)dV=limmax{ΔV}→0i=0∑nj=0∑mk=0∑tf(xi,yj,zk)dxidyjdzk

dx = input('dx = ');
dy = dx;dz = dx;
x = 0:dx:1;
for k = 1:length(x)x1 = (k-1)*dx;y = 0:dy:x1;for j = 1:length(y)y1 = (j-1)*dy;z1  =0:dz:x1.*y1;f = x1.*y1.^2.*z1.^3;s1(j) = trapz(f)*dz;ends2(k) = trapz(s1)*dy;
end
s = trapz(s2)*dx;
disp('值是:');
disp(s);

不知道为什么最后的结果和他算的不一样？？？

3.蒙特卡洛法

见博客通俗易懂的Monte Carlo积分方法（二）