题目六解答

问题1：非线性方程组求解

题目分析

这道题要求我们求方程组在(2,2)附近的根：

x² - y - 2 = 0
y² - 2x - 4 = 0

看到这个题目，我首先想到的是：这是一个非线性方程组，不能直接用线性代数的方法求解。我们需要用数值方法来逼近解。

解题思路

我选择使用牛顿法，因为：

牛顿法收敛速度快
对于这种光滑的非线性函数效果很好
我们有初始猜测值(2,2)，这给了我们一个好的起点

牛顿法的核心思想很简单：在初始点附近，我们用线性函数来近似非线性函数，然后逐步迭代逼近真实解。

代码实现

让我写一个牛顿法的函数：

function newton_method()
    # 从(2,2)开始
    x = 2.0
    y = 2.0
    
    # 定义我们的两个方程
    function f1(x, y)
        return x^2 - y - 2
    end
    
    function f2(x, y)
        return y^2 - 2*x - 4
    end
    
    # 计算雅可比矩阵 - 这是牛顿法的关键
    function jacobian(x, y)
        return [2*x -1; -2 2*y]
    end
    
    # 开始迭代
    for i in 1:10
        # 计算当前点的函数值
        f = [f1(x, y), f2(x, y)]
        # 计算雅可比矩阵
        J = jacobian(x, y)
        
        # 求解线性方程组 J * Δ = -f
        delta = J \ (-f)
        
        # 更新x和y
        x += delta[1]
        y += delta[2]
        
        # 如果变化很小，说明收敛了
        if abs(delta[1]) < 1e-10 && abs(delta[2]) < 1e-10
            break
        end
    end
    
    return x, y
end

结果

运行代码后，我得到了：

x = 2.21432
y = 2.903212

让我验证一下这个解是否正确：

x² - y - 2 = 0.0 ✓
y² - 2x - 4 = 0.0 ✓

解是正确的。

问题2：三重积分计算质量

题目分析

这道题比较复杂，要求我们计算一个复合零件的总质量。零件由三个曲面围成：

上部：抛物面 z = 4 - x² - y²
下部：圆锥面 z = √(x² + y²)
侧面：圆柱面 x² + y² = 1
密度函数：ρ(x, y, z) = 2 + z

解题思路

看到这个题目，我首先想到的是：这是一个三重积分问题。但是直接用直角坐标会很复杂，因为积分区域很特殊。

我注意到这个零件具有圆柱对称性（侧面是圆柱面），所以用柱坐标会大大简化计算。

坐标变换

在柱坐标下：

x = r·cos(θ)
y = r·sin(θ)
z = z

雅可比行列式是 r。

积分区域分析

让我分析一下积分区域：

r的范围：0到1（因为圆柱半径是1）
θ的范围：0到2π（完整的一圈）
z的范围：从圆锥面到抛物面，即从r到4-r²

代码实现

我写了两套代码：一套用解析积分，一套用数值积分来验证。

解析解

function calculate_mass_analytical()
    # 解析积分 ∫∫∫ (2 + z) r dr dθ dz
    
    # 先对z积分：∫(2r + rz) dz = 2r(z) + r(z²/2) | from r to 4-r²
    # 经过计算得到：16r - 6.5r³ - 2r² + 0.5r⁵
    
    # 然后对θ积分：2π * (16r - 6.5r³ - 2r² + 0.5r⁵)
    
    # 最后对r积分：∫(16r - 6.5r³ - 2r² + 0.5r⁵) dr from 0 to 1
    # 结果是：8 - 1.625 - 0.667 + 0.0833 = 5.7913
    
    mass_analytical = 2π * (8 - 1.625 - 2/3 + 1/12)
    return mass_analytical
end

数值解（用于验证）

function calculate_mass_numerical()
    # 用数值积分验证解析解
    n_r = 1000
    n_theta = 1000
    n_z = 1000
    
    mass = 0.0
    
    # 三重循环进行数值积分
    for i in 1:n_r
        r = (i - 0.5) / n_r  # 中点法则
        
        for j in 1:n_theta
            theta = 2π * (j - 0.5) / n_theta
            
            z_lower = r
            z_upper = 4 - r^2
            
            if z_upper > z_lower
                for k in 1:n_z
                    z = z_lower + (z_upper - z_lower) * (k - 0.5) / n_z
                    
                    # 密度函数乘以雅可比行列式
                    integrand = (2 + z) * r
                    
                    # 积分元素
                    dr = 1.0 / n_r
                    dtheta = 2π / n_theta
                    dz = (z_upper - z_lower) / n_z
                    
                    mass += integrand * dr * dtheta * dz
                end
            end
        end
    end
    
    return mass
end

结果

运行代码后，我得到了：

解析解：M = 36.390115 g
数值解：M = 36.390120 g

两个结果几乎完全一致，误差只有0.000005，这验证了我的解析解是正确的。

我还计算了体积和平均密度：

体积：V = 8.901180 cm³
平均密度：ρ_avg = 4.088235 g/cm³

验证过程

让我验证一下积分区域是否正确：

在r=0处：z从0到4 ✓
在r=1处：z从1到3 ✓
在r=0.5处：z从0.5到3.75 ✓

所有边界条件都满足，说明我的积分区域设置是正确的。

总结

问题1展示了如何用数值方法解决非线性方程组。牛顿法虽然简单，但非常有效，特别是当我们有好的初始猜测值时。

问题2是一个典型的三重积分问题，关键是要选择合适的坐标系。柱坐标让原本复杂的积分变得相对简单。我还用数值积分验证了解析解，这让我对结果更有信心。

下面是运行结果

下面是完整代码


println("题目六解答 ")
println("="^50)

# 问题1：求方程组在(x,y) = (2,2)附近的根
println("问题1：求方程组在(x,y) = (2,2)附近的根")
println("方程组：")
println("x² - y - 2 = 0")
println("y² - 2x - 4 = 0")
println()

# 使用牛顿法求解非线性方程组
function newton_method()
    # 初始猜测值
    x = 2.0
    y = 2.0
    
    # 定义函数
    function f1(x, y)
        return x^2 - y - 2
    end
    
    function f2(x, y)
        return y^2 - 2*x - 4
    end
    
    # 雅可比矩阵
    function jacobian(x, y)
        return [2*x -1; -2 2*y]
    end
    
    # 牛顿迭代
    for i in 1:10
        f = [f1(x, y), f2(x, y)]
        J = jacobian(x, y)
        
        # 求解线性方程组 J * Δ = -f
        delta = J \ (-f)
        
        x += delta[1]
        y += delta[2]
        
        # 检查收敛性
        if abs(delta[1]) < 1e-10 && abs(delta[2]) < 1e-10
            break
        end
    end
    
    return x, y
end

x_solution, y_solution = newton_method()

println("解：")
println("x = $(round(x_solution, digits=6))")
println("y = $(round(y_solution, digits=6))")
println()

# 验证解
println("验证解：")
println("x² - y - 2 = $(round(x_solution^2 - y_solution - 2, digits=10))")
println("y² - 2x - 4 = $(round(y_solution^2 - 2*x_solution - 4, digits=10))")

println()
println("="^50)

# 问题2：计算复合零件的总质量
println("问题2：计算复合零件的总质量")
println("零件边界：")
println("上部边界：抛物面 z = 4 - x² - y²")
println("下部边界：圆锥面 z = √(x² + y²)")
println("侧面边界：圆柱面 x² + y² = 1")
println("密度函数：ρ(x, y, z) = 2 + z (g/cm³)")
println()

# 使用柱坐标计算质量
# 在柱坐标下：
# x = r*cos(θ), y = r*sin(θ), z = z
# 密度函数：ρ(r, θ, z) = 2 + z
# 积分区域：0 ≤ r ≤ 1, 0 ≤ θ ≤ 2π, r ≤ z ≤ 4 - r²

function calculate_mass_analytical()
    # 解析积分
    # ∫∫∫ (2 + z) r dr dθ dz
    # = ∫∫∫ (2r + rz) dr dθ dz
    
    # 先对z积分：∫(2r + rz) dz = 2r(z) + r(z²/2) | from r to 4-r²
    # = 2r(4-r²-r) + r((4-r²)²/2 - r²/2)
    # = 2r(4-r²-r) + r((16-8r²+r⁴)/2 - r²/2)
    # = 2r(4-r²-r) + r(8-4r²+r⁴/2 - r²/2)
    # = 8r-2r³-2r² + 8r-4r³+r⁵/2 - r³/2
    # = 16r - 6.5r³ - 2r² + 0.5r⁵
    
    # 然后对θ积分：2π * (16r - 6.5r³ - 2r² + 0.5r⁵)
    
    # 最后对r积分：∫(16r - 6.5r³ - 2r² + 0.5r⁵) dr from 0 to 1
    # = [8r² - 1.625r⁴ - 0.667r³ + 0.0833r⁶] | from 0 to 1
    # = 8 - 1.625 - 0.667 + 0.0833 = 5.7913
    
    mass_analytical = 2π * (8 - 1.625 - 2/3 + 1/12)
    return mass_analytical
end

function calculate_mass_numerical()
    # 使用数值积分
    n_r = 1000
    n_theta = 1000
    n_z = 1000
    
    mass = 0.0
    
    for i in 1:n_r
        r = (i - 0.5) / n_r  # 中点法则
        
        for j in 1:n_theta
            theta = 2π * (j - 0.5) / n_theta
            
            z_lower = r
            z_upper = 4 - r^2
            
            if z_upper > z_lower
                for k in 1:n_z
                    z = z_lower + (z_upper - z_lower) * (k - 0.5) / n_z
                    
                    # 密度函数乘以雅可比行列式
                    integrand = (2 + z) * r
                    
                    # 积分元素
                    dr = 1.0 / n_r
                    dtheta = 2π / n_theta
                    dz = (z_upper - z_lower) / n_z
                    
                    mass += integrand * dr * dtheta * dz
                end
            end
        end
    end
    
    return mass
end

# 计算质量
mass_analytical = calculate_mass_analytical()
mass_numerical = calculate_mass_numerical()

println("零件的总质量：")
@printf("解析解：M = %.6f g\n", mass_analytical)
@printf("数值解：M = %.6f g\n", mass_numerical)
println()

# 计算体积
function calculate_volume()
    # 体积 = ∫∫∫ r dr dθ dz
    volume = 0.0
    
    n_r = 1000
    n_theta = 1000
    n_z = 1000
    
    for i in 1:n_r
        r = (i - 0.5) / n_r
        
        for j in 1:n_theta
            theta = 2π * (j - 0.5) / n_theta
            
            z_lower = r
            z_upper = 4 - r^2
            
            if z_upper > z_lower
                for k in 1:n_z
                    z = z_lower + (z_upper - z_lower) * (k - 0.5) / n_z
                    
                    # 雅可比行列式
                    integrand = r
                    
                    # 积分元素
                    dr = 1.0 / n_r
                    dtheta = 2π / n_theta
                    dz = (z_upper - z_lower) / n_z
                    
                    volume += integrand * dr * dtheta * dz
                end
            end
        end
    end
    
    return volume
end

volume = calculate_volume()
println("零件体积：")
@printf("V = %.6f cm³\n", volume)

# 平均密度
avg_density = mass_analytical / volume
println("平均密度：")
@printf("ρ_avg = %.6f g/cm³\n", avg_density)

println()
println("积分区域验证：")
println("在r=0处：z从0到4")
println("在r=0.5处：z从0.5到3.75")
println("在r=1处：z从1到3")
println("在r=1.5处：z从1.5到1.75")
println("在r=2处：z从2到0（超出边界）")

# 验证解析解的正确性
println()
println("解析解验证：")
println("对z积分：∫(2+z)dz = 2z + z²/2 | from r to 4-r²")
println("= 2(4-r²) + (4-r²)²/2 - 2r - r²/2")
println("= 8-2r² + (16-8r²+r⁴)/2 - 2r - r²/2")
println("= 8-2r² + 8-4r²+r⁴/2 - 2r - r²/2")
println("= 16 - 6.5r² - 2r + 0.5r⁴")

题目六解答

问题1：非线性方程组求解

题目分析

解题思路

代码实现

结果

问题2：三重积分计算质量

题目分析

解题思路

坐标变换

积分区域分析

代码实现

解析解

数值解（用于验证）

结果

验证过程

总结

全部回答