MATLAB 数学实验

2 微分方程^[1]

针对不同的阻尼系数用 ODE 求数值解

自然频率 \(\omega=10\)，初始条件 \(x(0)=1,x'(0)=0\)。求解

\[ \begin{gather} \ddot x+20c\dot x+100x=0. \end{gather} \]

function dx = odefun4(t, x)
    global c;
    dx = zeros(2, 1);
    dx(1) = x(2);
    dx(2) = -20 * c * x(2) - 100 * x(1);
end

% ode45(odefun, tspan, y0)
% odefun: 函数句柄
% tspan: 自变量区间
% y0: 初值
vc = [0.2 0.4 0.6 0.8 1];
tspan = linspace(0, 4, 100);
x0 = [1 0];
global c;
hold on;

for i = 1:length(vc)
    c = vc(i);
    [t, x] = ode45('odefun4', tspan, x0);
    text(t(10), x(10, 1), ['\leftarrow c=', num2str(c)]);
    plot(t, x(:, 1)); % 解的函数图
end

可使用 plot(t, x(:, 2)) 绘制解的导数图。

弹簧振动动图

vc = [0.2 0.4 0.6 0.8 1];
tspan = linspace(0, 4, 200);
x0 = [1 0];
global c;
c = 0.1;
hold on;

[t, xt] = ode45('odefun4', tspan, x0);
animinit('onecart1 Animation')
axis([-2 6 -10 10]);
u = 2;
xy = [0 0 0 0 u u u + 1 u + 1 u u; -1.2 0 1.2 0 0 1.2 1.2 -1.2 -1.2 0];
x = xy(1, :);
y = xy(2, :);
plot([-10 20], [-1.4 -1.4], 'b-', 'LineWidth', 2);
hndl = plot(x, y, 'b-', 'EraseMode', 'XOR', 'LineWidth', 2);
set(gca, 'UserData', hndl);

for i = 1:length(t)
    u = 2 + 5 * xt(i);
    x = [0 0 0 0 u u u + 1 u + 1 u u];
    hndl = get(gca, 'UserData');
    set(hndl, 'XData', x);
    pause(0.02);
    drawnow;
end

---

分形实例

Koch 曲线

% 参数设置
k = 5; % 迭代次数
n = 1; % 线段数量
theta = pi / 3; % 旋转角度

% 初始线段端点坐标
p = [0 0; 10 0];

% 旋转矩阵
A = [cos(theta) -sin(theta); sin(theta) cos(theta)];

% 存放新的线段坐标
r = zeros(4 * n, 2);

% 迭代处理
for iter = 1:k
    j = 0;

    for i = 1:n
        % 当前线段的起点和终点
        q1 = p(i, :);
        q2 = p(i + 1, :);

        % 计算新点坐标
        d = (q2 - q1) / 3;
        new_points = [q1; q1 + d; q1 + d + d * A'; q1 + 2 * d];

        % 存入r
        r(j + 1:j + 4, :) = new_points;
        j = j + 4;
    end

    % 更新线段数量和端点坐标
    n = 4 * n;
    p = [r; q2]; % 重新装载本次迭代的所有点
end

plot(p(:, 1), p(:, 2));
fill(p(:, 1), p(:, 2), 'k'); % 填充颜色
set(findobj(gcf, 'type', 'patch'), 'edgecolor', 'none'); % 去掉边框
axis off;
axis equal;

Sierpinski 地毯

x = 0; % 初始正方形左下角顶点坐标
y = 0;
d = 1; % 正方形边长
n = 3; % 迭代次数

% 迭代处理
for p = 1:n
    % 迭代后所有小正方形左下角顶点坐标
    a1 = [];
    b1 = [];

    for q = 1:length(x) % 对每个小正方形进行迭代
        % 新的小正方形左下角顶点坐标
        x1 = x(q) + [0 d / 3 2 * d / 3 0 2 * d / 3 0 d / 3 2 * d / 3];
        y1 = y(q) + [0 0 0 d / 3 d / 3 2 * d / 3 2 * d / 3 2 * d / 3];
        a1 = [a1 x1];
        b1 = [b1 y1];
    end

    d = d / 3;
    x = a1;
    y = b1;
end

hold on;

for q = 1:length(x)
    fill(x(q) + [0 d d 0 0], y(q) + [0 0 d d 0], 'b');
end

set(findobj(gcf, 'type', 'patch'), 'edgecolor', 'none'); % 去掉边框
axis off;
axis equal;

分形树木

function [z, A] = frat3(z, A, n)
    % z: 起始点
    % A: 起始角度
    % n: 迭代次数
    N = 6; % 最大迭代次数
    s = 0.7; % 缩放因子
    Br = [pi / 10; -pi / 10]; % 旋转角度

    if n > N
        return
    end

    hold on;
    k = 1;

    if n == 1 % 第一次迭代，画主干
        plot([z z + exp(i * A)], 'linewidth', 3 / N * (N - n + 1));
        pause(0.05);
        [z, A] = frat3(z + exp(i * A), A, n + 1); % 画分支
    else

        for k = 1:2 % k=1,2分别画两个分支
            A = A + Br(k);
            leng = s ^ (n - 1);
            plot([z z + leng * exp(i * A)], 'linewidth', 3 / N * (N - n + 1));
            pause(0.05);
            z = z + leng * exp(i * A);
            [z, A] = frat3(z, A, n + 1); % 第n+1次迭代
            A = A + pi;
            z = z + leng * exp(i * A);
            A = A + pi - Br(k);
            k = k + 1;
        end

    end

end

使用 frat3(0, pi/2, 1) 调用。

参考文献

1. 黄平，刘小兰，温旭辉. 数学实验典型案例. 高等教育出版社, 2015. 2019 年 11 月第 4 次印刷. ↩︎