十、例題

【例3】求由彈簧連接的 $100$ 個質點的位移 $u(1),u(2),...,u(100)$, 彈性系數均為 $c=1$, 每個質點受到的外力均為 $f(i)=0.01$. 畫出兩端固定和固定-自由這兩種情形 u 的圖形。
解：

% 參數設置
n = 100;            % 質點的數量
c = 1;              % 彈性系數
f = 0.01*ones(n,1); % 每個質點受到的外力% 1. 兩端固定邊值條件
% 使用diag函數構造矩陣A
A0 = diag(ones(n+1, 1), 0)-diag(ones(n,1), -1);        % A0是一個101*100的矩陣，這里創建一個101*101的矩陣
A0 = A0(:, 1:n);                                       % 截取前100列即可得到矩陣A0% 構造剛度矩陣K
K_fixed = A0'*A0;                   % 由于彈性系數均為1，所以此處不含有矩陣C% 求解位移u
u_fixed = K_fixed \ f;% 繪制兩端固定邊值條件的位移圖
figure(1);
plot(u_fixed, '+');
xlabel('mass number');           % 質點編號
ylabel('movement');              % 位移
title('兩端固定邊值條件的位移'); 
grid on;
print -dpng fixed_fixed_displacement.png;% 2. 固定-自由邊值條件
% 構造矩陣 A
A1 = A0(1:n, :);        % A1是一個方陣，取A0的前100行即可% 構造剛度矩陣K
K_free = A1'*A1;% 求解位移
u_free = K_free \ f;% 繪制固定-自由邊值條件的位移圖
figure(2);
plot(u_free, '+');
xlabel('mass number');
ylabel('movement');
title('固定-自由邊值條件的位移');
grid on;
print -dpng fixed_free_displacement.png

運行結果：

【例4】化學工程中通常用一階導數 $\frac{\text{d}u}{\text{d}x}$ 表示流體流速，用二階導數 $\frac{{\text{d}}^{2}u}{\text{d}{x}^2}$ 描述擴散速度。將 $\frac{\text{d}u}{\text{d}x}$ 分別換成向前、中心差分和向后差分，取 $\Delta x=\frac{1}{8}$. 畫出這三種離散方法求出的下面方程數值解的圖形：
$$?\frac{{\text{d}}^{2}u}{\text{d}{x}^2}+10\frac{\text{d}u}{\text{d}x}=1,u(0)=u(1)=0.$$
這種 對流-擴散方程(convection-diffsion equation) 無處不在，這個方程轉化為描述期權定價問題的布萊克-斯科爾斯（Black-Scholes)方程。
解：

E = diag(ones(6,1),1);      % 創建一個對角矩陣，用來輔助計算 K 和 D
K = 64 * (2*eye(7)-E-E');   % 計算二階差分矩陣
D = 80 * (E-eye(7));        % 向前一階差分矩陣
u = (K+D)\ones(7,1);        % 向前差分求解x = linspace(0, 1, 9);      % 區間[0,1]之間等間距9個點，則每段長度1/8  % 畫出圖形
figure;
plot(x, [0; u; 0], 'o-', 'DisplayName', 'Forward Difference');  % 繪制向前差分解的圖像并標上圖例, [0; u; 0] 是拼接列向量
hold on;
u = (K-D')\ones(7,1);           % 向后差分
plot(x, [0; u; 0], 'd-', 'DisplayName', 'Backward Difference'); % 繪制向后差分解的圖像并標上圖例
u = (K+D/2-D'/2)\ones(7,1);     % 中心差分
plot(x, [0; u; 0], 's-', 'DisplayName', 'Centred Difference');  % 繪制中心差分解的圖像并標上圖例
xlabel('x'); ylabel('u(x)');    % 坐標軸表示
title('Numrerical Solutions of Convection-Diffusion Equation'); % 圖像標題
legend('Location', 'northwest');                                % 圖例在左上角
grid on;% 精確解計算：u(x) = 1/(10(e^{10}-1))(1-e^{10x})+x/10
e10 = exp(10);
A = 1/(10*(e10-1));
B = -1/(10*(e10-1));
u_exact = A + B * exp(10*x) + x/10;
plot(x, u_exact, 'k-', 'DisplayName', 'Exact Solution');         % 畫出第2~8個點的精確解圖像
hold off;

運行結果：

通常情況下，中心差分是最優解，它最接近精確解。