方程求解

多项式求根

roots 函数

使用 roots 函数求解多项式方程的根。多项式由系数向量表示，从最高次到常数项：

# 多项式：p(x) = a_n*x^n + ... + a_1*x + a_0
# 系数向量：[a_n, ..., a_1, a_0]

r = roots(coeffs)

一次方程

# 求解：2x + 3 = 0
coeffs = [2, 3]
r = roots(coeffs)
disp("解：x = " + num2str(r))  # -1.5

二次方程

# 求解：x^2 - 5x + 6 = 0
# (x-2)(x-3) = 0，根为 2 和 3
coeffs = [1, -5, 6]
r = roots(coeffs)
disp("根：")
disp(r)  # [3; 2]

# 求解：x^2 + 2x + 5 = 0（有复数根）
coeffs = [1, 2, 5]
r = roots(coeffs)
disp("复数根：")
disp(r)  # [-1+2i; -1-2i]

高次方程

# 求解：x^3 - 6x^2 + 11x - 6 = 0
# (x-1)(x-2)(x-3) = 0
coeffs = [1, -6, 11, -6]
r = roots(coeffs)
disp("三次方程的根：")
disp(r)  # [3; 2; 1]

# 求解：x^4 - 1 = 0
coeffs = [1, 0, 0, 0, -1]
r = roots(coeffs)
disp("四次方程的根：")
disp(r)  # [1; -1; i; -i]

从根构造多项式

# poly 函数：从根生成多项式系数
roots_vec = [1, 2, 3]
coeffs = poly(roots_vec)
disp("多项式系数：")
disp(coeffs)  # [1 -6 11 -6]

# 验证
r = roots(coeffs)
disp("验证根：")
disp(r)

线性方程组

使用左除运算符 \

求解线性方程组 Ax = b，使用矩阵左除 \：

# Ax = b
x = A \ b

基本示例

# 求解方程组：
#   2x + y = 5
#   x - y = 1
# 即：x = 2, y = 1

A = [2 1; 1 -1]
b = [5; 1]

x = A \ b
disp("解：")
disp(x)  # [2; 1]

# 验证
result = A * x
disp("验证 Ax = b：")
disp(result)  # [5; 1]

三元方程组

# 求解：
#   x + 2y + z = 6
#   2x + y - z = 1
#   x - y + 2z = 5

A = [1 2 1; 2 1 -1; 1 -1 2]
b = [6; 1; 5]

x = A \ b
disp("解向量：")
disp(x)

# 提取各个变量
x1 = x(1)
x2 = x(2)
x3 = x(3)
disp("x = " + num2str(x1))
disp("y = " + num2str(x2))
disp("z = " + num2str(x3))

欠定和超定系统

# 超定系统（方程多于未知数）- 最小二乘解
A = [1 1; 1 2; 1 3]  # 3个方程，2个未知数
b = [2; 3; 4]

x = A \ b  # 自动计算最小二乘解
disp("最小二乘解：")
disp(x)

# 计算残差
residual = b - A * x
disp("残差：")
disp(residual)

齐次方程组

# 求解 Ax = 0 的非零解
A = [1 2 3; 4 5 6; 7 8 9]  # 秩不满

# 检查秩
r = rank(A)
disp("矩阵的秩：" + num2str(r))

# 对于齐次方程，可以通过求零空间获得解
# 使用 SVD 分解等方法（如果支持）

非线性方程求解

使用 solve 函数

TianYuan 提供 solve 函数求解非线性方程：

# solve(f, x0, tolerance, maxIter)
# f: 函数
# x0: 初始猜测
# tolerance: 容差
# maxIter: 最大迭代次数

root = solve(@myFunc, x0)

示例：求平方根

# 定义方程 f(x) = x^2 - 2 = 0
function y = f(x)
    y = x^2 - 2
end

# 求解（从初始猜测 1.0 开始）
root = solve(@f, 1.0)
disp("sqrt(2) ≈ " + num2str(root))  # 1.414...

示例：超越方程

# 求解：x = cos(x)
function y = f(x)
    y = x - cos(x)
end

root = solve(@f, 0.5)
disp("x = cos(x) 的解：" + num2str(root))  # 0.739...

手动实现牛顿法

# 牛顿法求根：x_{n+1} = x_n - f(x_n)/f'(x_n)
function root = newtonMethod(f, df, x0, tol, maxIter)
    x = x0
    
    for iter = 1:maxIter
        fx = f(x)
        
        if abs(fx) < tol
            root = x
            return
        end
        
        dfx = df(x)
        x = x - fx / dfx
    end
    
    root = x
    disp("警告：达到最大迭代次数")
end

# 示例：求 sqrt(2)
function y = f(x)
    y = x^2 - 2
end

function y = df(x)
    y = 2*x
end

root = newtonMethod(@f, @df, 1.0, 1e-10, 100)
disp("sqrt(2) = " + num2str(root))

常微分方程（ODE）

ode45 函数

使用 Dormand-Prince RK45 自适应步长方法求解常微分方程初值问题：

# ode45(f, tspan, y0)
# f: 函数句柄 @f，其中 f(t, y) 返回 dydt
# tspan: 时间区间 [t0, tf]
# y0: 初始条件（标量或列向量）
# 返回: t 为列向量，y 为矩阵（每行是一个时刻的状态）

[t, y] = ode45(@odefun, tspan, y0)

说明：默认相对容差 1e-3，绝对容差 1e-6。对于刚性方程组可能出现步长过小的报错， 此类方程建议手动实现隐式方法（如 Backward Euler）。

关键概念：返回值结构

ode45 返回两个输出：t 是列向量（各时刻），y 是矩阵（每行对应一个时刻，每列对应一个状态分量）：

# 一阶 ODE — 标量状态
function dydt = decay(t, y)
    dydt = -2 * y        # dy/dt = -2y
end

[t, y] = ode45(@decay, [0, 3], 1)

# t: 列向量，形如 [0; 0.05; 0.12; ...]
# y: 同长度的列向量，y(k,1) 是 t(k) 时刻的状态值

disp("步数 = " + num2str(length(t)))     # 自适应步数（约数十步）
disp("y(0) = " + num2str(y(1, 1)))       # 1.0（初始值）
disp("y(3) ≈ " + num2str(y(length(t), 1))) # ≈ 0.00248（e^{-6}）

关键概念：向量状态（方程组）

对于 n 维方程组，y0 是 n×1 列向量，ODE 函数返回同形状的列向量，y 结果矩阵有 n 列：

# 二阶 ODE 转为一阶方程组：d²x/dt² + x = 0
# 状态 y = [x; x']，方程组 y' = [y(2); -y(1)]
function dydt = harmonic(t, y)
    dydt = [y(2); -y(1)]   # 列向量返回
end

y0 = [1; 0]                # 初始位移 1，初始速度 0
[t, y] = ode45(@harmonic, [0, 6.2832], y0)

# y 是 N×2 矩阵：
#   第 1 列 = 位移 x(t)
#   第 2 列 = 速度 x'(t)
x = y(:, 1)   # 提取所有时刻的位移
v = y(:, 2)   # 提取所有时刻的速度

plot(t, x)
xlabel('t'); ylabel('x(t)')
title('简谐振动（ω=1）')

示例：简单的一阶ODE

# 求解：dy/dt = -2y, y(0) = 1
# 解析解：y(t) = e^(-2t)

function dydt = odefun(t, y)
    dydt = -2 * y
end

# 在 [0, 2] 区间求解
tspan = [0, 2]
y0 = 1

[t, y] = ode45(@odefun, tspan, y0)

# 绘制结果
plot(t, y)
xlabel('t')
ylabel('y')
title('dy/dt = -2y 的解')
grid on

# 与解析解比较
y_exact = exp(-2 * t)
hold on
plot(t, y_exact, '--')
legend('数值解', '解析解')
hold off

示例：指数增长模型

# 人口增长模型：dP/dt = rP
# P(0) = P0, r = 0.1

function dPdt = population(t, P)
    r = 0.1
    dPdt = r * P
end

tspan = [0, 50]
P0 = 100

[t, P] = ode45(@population, tspan, P0)

plot(t, P)
xlabel('时间')
ylabel('人口')
title('指数增长模型')
grid on

示例：简谐振动

# 二阶ODE：d²x/dt² + ω²x = 0
# 转换为一阶系统：
# dx/dt = v
# dv/dt = -ω²x

function dydt = harmonic(t, y)
    omega = 2 * pi  # 角频率
    
    x = y(1)
    v = y(2)
    
    dxdt = v
    dvdt = -omega^2 * x
    
    dydt = [dxdt; dvdt]
end

# 初始条件：x(0) = 1, v(0) = 0
tspan = [0, 2]
y0 = [1; 0]

[t, y] = ode45(@harmonic, tspan, y0)

# 提取位置和速度
x = y(:, 1)
v = y(:, 2)

# 绘制位置
plot(t, x)
xlabel('时间')
ylabel('位移')
title('简谐振动')
grid on

示例：洛伦兹吸引子

# 洛伦兹方程（混沌系统）：
# dx/dt = σ(y - x)
# dy/dt = x(ρ - z) - y
# dz/dt = xy - βz

function dydt = lorenz(t, y)
    sigma = 10
    rho = 28
    beta = 8/3
    
    x = y(1)
    y_val = y(2)
    z = y(3)
    
    dxdt = sigma * (y_val - x)
    dydt_val = x * (rho - z) - y_val
    dzdt = x * y_val - beta * z
    
    dydt = [dxdt; dydt_val; dzdt]
end

# 初始条件
tspan = [0, 50]
y0 = [1; 1; 1]

[t, y] = ode45(@lorenz, tspan, y0)

# 绘制三维轨迹
plot3(y(:,1), y(:,2), y(:,3))
xlabel('X')
ylabel('Y')
zlabel('Z')
title('洛伦兹吸引子')
grid on

优化问题

黄金分割搜索（一维优化）

function xmin = goldenSearch(f, a, b, tol)
    # 黄金分割法找函数最小值
    phi = (1 + sqrt(5)) / 2  # 黄金比例
    
    while abs(b - a) > tol
        x1 = b - (b - a) / phi
        x2 = a + (b - a) / phi
        
        if f(x1) < f(x2)
            b = x2
        else
            a = x1
        end
    end
    
    xmin = (a + b) / 2
end

# 示例：找 f(x) = (x-2)^2 的最小值
function y = parabola(x)
    y = (x - 2)^2
end

xmin = goldenSearch(@parabola, 0, 5, 1e-6)
disp("最小值点：x = " + num2str(xmin))  # 2
disp("最小值：f(x) = " + num2str(parabola(xmin)))  # 0

梯度下降（多维优化）

function xmin = gradientDescent(f, grad, x0, alpha, tol, maxIter)
    # 梯度下降法
    # f: 目标函数
    # grad: 梯度函数
    # x0: 初始点
    # alpha: 学习率
    # tol: 容差
    # maxIter: 最大迭代次数
    
    x = x0
    
    for iter = 1:maxIter
        g = grad(x)
        
        if norm(g) < tol
            xmin = x
            return
        end
        
        x = x - alpha * g
    end
    
    xmin = x
end

# 示例：最小化 f(x,y) = x^2 + y^2
function z = f(x)
    z = x(1)^2 + x(2)^2
end

function g = grad(x)
    g = [2*x(1); 2*x(2)]
end

x0 = [5; 5]
xmin = gradientDescent(@f, @grad, x0, 0.1, 1e-6, 1000)
disp("最小值点：")
disp(xmin)  # [0; 0]

综合示例

示例1：求方程组的交点

# 求两个圆的交点
# 圆1: x^2 + y^2 = 1
# 圆2: (x-1)^2 + y^2 = 1

# 从第二个方程减去第一个方程：
# -2x + 1 = 0 => x = 0.5
# 代入第一个方程：0.25 + y^2 = 1 => y = ±sqrt(0.75)

x = 0.5
y = sqrt(0.75)

disp("交点1：(" + num2str(x) + ", " + num2str(y) + ")")
disp("交点2：(" + num2str(x) + ", " + num2str(-y) + ")")

示例2：弹道轨迹

# 在重力作用下的抛体运动
# d²x/dt² = 0
# d²y/dt² = -g

function dydt = projectile(t, y)
    g = 9.8  # 重力加速度
    
    # y = [x, vx, y_pos, vy]
    x = y(1)
    vx = y(2)
    y_pos = y(3)
    vy = y(4)
    
    dxdt = vx
    dvxdt = 0
    dydt_pos = vy
    dvydt = -g
    
    dydt = [dxdt; dvxdt; dydt_pos; dvydt]
end

# 初始条件：v0 = 20 m/s, θ = 45°
v0 = 20
theta = 45 * pi / 180

vx0 = v0 * cos(theta)
vy0 = v0 * sin(theta)

tspan = [0, 3]
y0 = [0; vx0; 0; vy0]

[t, y] = ode45(@projectile, tspan, y0)

# 提取轨迹
x = y(:, 1)
y_pos = y(:, 3)

# 只绘制在地面以上的部分
idx = y_pos >= 0
plot(x(idx), y_pos(idx))
xlabel('水平距离 (m)')
ylabel('高度 (m)')
title('抛体运动轨迹')
grid on

示例3：化学反应动力学

# 一级反应：A -> B
# d[A]/dt = -k[A]
# d[B]/dt = k[A]

function dcdt = reaction(t, c)
    k = 0.5  # 反应速率常数
    
    A = c(1)
    B = c(2)
    
    dAdt = -k * A
    dBdt = k * A
    
    dcdt = [dAdt; dBdt]
end

# 初始浓度：[A]0 = 1, [B]0 = 0
tspan = [0, 10]
c0 = [1; 0]

[t, c] = ode45(@reaction, tspan, c0)

# 绘制浓度随时间的变化
plot(t, c(:,1), 'b-', t, c(:,2), 'r-')
xlabel('时间')
ylabel('浓度')
legend('[A]', '[B]')
title('化学反应动力学')
grid on
下一章：图形绘制 →