美賽整理之Matlab的工程數學計算學習筆記（高等數學）

1.極限的定義和判別：

​ 由極限的定義可知：對於給定的$f(x)$,$\forall ϵ>0,\exists \delta >0，$使得：當$|x-x_0|<\delta$時，有$|f(x)-A|<ϵ$成立，那麼我們稱：$li{m}_{x\to x_0}f(x)=A$。

​ 那麼我們的關鍵就是要給出$\forall ϵ>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 == 1
    epsilon = input('請輸入一個小鼠:\n');
    while abs(A - eval(fxc))>epsilon 
        delta = delta/2;
       
       if p == 1 
           x = x_c - delta;
       else
           x = x_c + delta;
       end
       
        if abs(delta)<eps
            disp('沒有找到這樣的delta？');
            n = 0;
            break
        end
    end
    if n == 0
        if p == 1
           disp('左極限不正確。');
        else
           disp('右極限不正確');
        end
        break
    end
    if n==1
        if p == 1
            disp('左極限可能正確');
        else
            disp('右極限可能正確');
        end
        disp('delta的值是：');
        disp(delta);
    end
    flag = input('還要重新試一下嗎？1：重新試一下/0：不試了\n');
    if flag == 1
       p = input('請告訴我這次要判斷左極限還是右極限？1：左極限，0：右極限\n');
    end
end

2.繪製特殊曲面

若二次曲面的方程如下：
$$\frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2}=d$$
根據$a,b,c$的不同繪製出不同的圖形。

在給定$x,y$要求求出$z$時，如果有開方的運算，會出現虛數，而且對於實數也會有兩個符號不同的解。為了使得虛數不會出現在繪圖中，我們把虛數換成·非數($NaN$)

我們規定：如果一個z值只有只有這一點時虛數，我們就讓它為0。如果這個數四周都是虛數，我們就讓它時$NaN$

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-1
    for j = 2:N-1
        if(imag(z(i,j))~=0)
            z1(i,j) = 0;
        end
        if all(imag(z([i-1:i+1],[j-1:j+1]))~=0)
            z1(i,j) = NaN;
        end
    end
end
surf(x,y,z1);
hold on
if u == 1
    z2 = - 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.求兩個空間曲面的交線

求以下兩個方程的交線：
$$\begin{gathered} z=x^2-2y^2 \\ z=2x-3y \end{gathered}$$

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.m$函式。再呼叫Matlab的函式：$s=quad{(}^{\prime }exampl{e}^{\prime },a,b)$

5.多重積分的計算

1.截面法：

∫ ∫ D ( x 2 + y 2 ) 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:lx
    x1 = (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/3$。

2.定義法

計算：
$$\int \int {\int }_{\Omega }x{y}^2{z}^{3}dxdydz$$
其中：$\Omega$是由$x=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 ) d V = l i m m a x { Δ V } → 0 ∑ i = 0 n ∑ j = 0 m ∑ k = 0 t f ( x i , y j , z k ) d x i d y j d z k \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}​→0​i=0∑n​j=0∑m​k=0∑t​f(xi​,yj​,zk​)dxi​dyj​dzk​

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;
    end
    s2(k) = trapz(s1)*dy;
end
s = trapz(s2)*dx;
disp('值是:');
disp(s);

不知道為什麼最後的結果和他算的不一樣？？？

3.蒙特卡洛法

見部落格通俗易懂的Monte Carlo積分方法（二）