22、数值计算与编程实践：调试、求解与模拟

数值计算与编程实践：调试、求解与模拟

1. 调试技巧与函数修正

在编程中，调试是解决代码问题的关键步骤。以绘制两条曲线之间区域的函数 shadecurves 为例，我们需要对其进行修正以确保正常工作。

1.1 修正思路

原代码可能存在问题，通过对 x 坐标的处理进行调整，即反转 xvals 得到 slavx ，从而修正函数。

1.2 修正后的代码

function shadecurves(f, g, a, b)
%SHADECURVES Draws the region between two curves
% SHADECURVES(f, g, a, b) takes strings or expressions f
% and g, interprets them as functions, plots them between
% x = a and x = b, and shades the region in between.
% Example: shadecurves('sin(x)', '-sin(x)', 0, pi)
ffun = inline(vectorize(f)); gfun = inline(vectorize(g));
xvals = a:(b - a)/50:b; slavx = b:(a - b)/50:a;
patch([xvals,slavx], [ffun(xvals),gfun(slavx)], [.2,0,.8])

修正后的函数能够正确绘制两条曲线之间的区域，可对其他示例进行测试。

2. 实践集 A 解答

2.1 基本运算

• 减法运算 ： 1111 - 345 的结果为 766 。
• 指数与乘法运算 ：

format long; [exp(14), 382801*pi]

结果显示第二个数更大。
- 数值近似比较 ：

[2709/1024, 10583/4000, 2024/765, sqrt(7)]

其中 2024/765 是最佳近似值。

2.2 函数计算

函数	代码	结果
双曲余弦函数	cosh(0.1)	1.00500416805580
自然对数函数	log(2)	0.69314718055995
反正切函数	atan(1/2)	0.46364760900081

2.3 方程求解

• 三元一次方程组求解 ：

[x, y, z] = solve('3*x + 4*y + 5*z = 2', '2*x - 3*y + 7*z = -1', 'x - 6*y + z = 3', 'x', 'y', 'z')

求解结果可通过矩阵乘法进行验证。
- 另一个三元一次方程组求解 ：

[x, y, z] = solve('3*x - 9*y + 8*z = 2', '2*x - 3*y + 7*z = -1', 'x - 6*y + z = 3', 'x', 'y', 'z')

得到一个依赖于变量 y 的单参数解族，同样可通过矩阵乘法验证。

2.4 因式分解与化简

• 因式分解 ：

syms x y; factor(x^4 - y^4)

结果为 (x - y)*(x + y)*(x^2 + y^2) 。
- 表达式化简 ：

simplify(1/(1 + 1/(1 + 1/x)))
simplify(cos(x)^2 - sin(x)^2)

可使用 simple 函数进一步化简。

2.5 幂运算与符号计算

3^301
sym('3')^301
sym('3^301')

注意 sym 函数本身不会进行求值。

2.6 方程求解与绘图

• 线性方程求解 ：

solve('8*x + 3 = 0', 'x')
vpa(ans, 15)

• 三次方程求解 ：

syms p q; solve('x^3 + p*x + q = 0', 'x')

• 函数绘图与零点求解 ：

ezplot('exp(x)'); hold on; ezplot('8*x - 4'); hold off
fzero(inline('exp(x) - 8*x + 4'), 1)
fzero(inline('exp(x) - 8*x + 4'), 3)

2.7 函数绘图

• 多项式函数绘图 ：

ezplot('x^3 - x', [-4 4])

• 三角函数绘图 ：

ezplot('sin(1/x^2)', [-2 2])

由于函数在 x = 0 附近振荡剧烈，绘图可能不完整，可细化网格进行尝试。
- 其他函数绘图 ：

ezplot('tan(x/2)', [-pi pi]); axis([-pi, pi, -10, 10])
X = -2:0.05:2;
plot(X, exp(-X.^2), X, X.^4 - X.^2)

2.8 函数交点求解

通过绘图和 fzero 函数求解 2^x 和 x^4 的交点：

ezplot('x^4'); hold on; ezplot('2^x'); hold off
X = -2:0.1:2; plot(X, 2.^X); hold on; plot(X, X.^4, '--'); hold off
X = 6:0.1:20; plot(X, 2.^X); hold on; plot(X, X.^4, '--'); hold off
r1 = fzero(inline('2^x - x^4'), -0.9)
r2 = fzero(inline('2^x - x^4'), 1.2)
r3 = fzero(inline('2^x - x^4'), 16)

使用 solve 函数可得到包含复数解的符号解。

3. 实践集 B 解答

3.1 等高线绘图

graph LR
    A[定义网格] --> B[绘制等高线]
    B --> C[调整参数]
    C --> D[绘制特定等高线]

• 不同范围的等高线绘图 ：

[X, Y] = meshgrid(-1:0.1:1, -1:0.1:1); contour(X, Y, 3*Y + Y.^3 - X.^3, 'k')
[X, Y] = meshgrid(-10:0.1:10, -10:0.1:10); contour(X, Y, 3*Y + Y.^3 - X.^3, 'k')
[X, Y] = meshgrid(-10:0.1:10, -10:0.1:10); contour(X, Y, 3*Y + Y.^3 - X.^3, 30, 'k')

• 特定等高线绘图 ：

[X, Y] = meshgrid(-5:0.1:5, -5:0.1:5); contour(X, Y, 3.*Y + Y.^3 - X.^3, [5 5], 'k')

• 通过特定点的等高线绘图 ：

[X, Y] = meshgrid(0:0.1:2, 0:0.1:2); contour(X, Y, Y.*log(X) + X.*log(Y), [0 0], 'k')

3.2 导数计算

函数	代码	导数结果
多项式函数	diff(6*x^3 - 5*x^2 + 2*x - 3, x)	18*x^2 - 10*x + 2
分式函数	diff((2*x - 1)/(x^2 + 1), x)	可化简
三角函数	diff(sin(3*x^2 + 2), x)	6*cos(3*x^2 + 2)*x
反三角函数	diff(asin(2*x + 3), x)	1/(-2 - x^2 - 3*x)^(1/2)
根式函数	diff(sqrt(1 + x^4), x)	2/(1 + x^4)^(1/2)*x^3
幂函数	diff(x^r, x)	x^r*r/x
反正切函数	diff(atan(x^2 + 1), x)	2*x/(1 + (x^2 + 1)^2)

3.3 积分计算

积分类型	代码	结果
定积分	int(cos(x), x, 0, pi/2)	1
不定积分	int(x*sin(x^2), x)	-1/2*cos(x^2)
换元积分	int(sin(3*x)*sqrt(1 - cos(3*x)), x)	2/9*(1 - cos(3*x))^(3/2)
根式积分	int(x^2*sqrt(x + 4), x)	可化简
广义积分	int(exp(-x^2), x, -Inf, Inf)	pi^(1/2)

3.4 数值积分

int(exp(sin(x)), x, 0, pi)
format long; quadl('exp(sin(x))', 0, pi)
quadl('sqrt(x.^3 + 1)', 0, 1)
quadl('exp(-x.^2)', -35, 35)

可与精确积分结果进行比较。

3.5 极限计算

limit(sin(x)/x, x, 0)
limit((1 + cos(x))/(x + pi), x, -pi)
limit(x^2*exp(-x), x, Inf)
limit(1/(x - 1), x, 1, 'left')
limit(sin(1/x), x, 0, 'right')

部分极限可通过绘图辅助理解。

3.6 级数求和

graph LR
    A[符号定义] --> B[级数求和]
    B --> C[化简或处理错误]

• 平方级数求和 ：

syms k n r x z
symsum(k^2, k, 0, n)
simplify(ans)

• 等比级数求和 ：

symsum(r^k, k, 0, n)

• 指数级数求和 ：

symsum(x^k/factorial(k), k, 0, Inf)

处理错误可使用 gamma 函数或 maple 函数。
- 其他级数求和 ：

symsum(1/(z - k)^2, k, -Inf, Inf)

3.7 泰勒展开

taylor(exp(x), 7, 0)
taylor(sin(x), 5, 0)
taylor(sin(x), 6, 0)
taylor(sin(x), 6, 2)
taylor(tan(x), 7, 0)
taylor(log(x), 5, 1)
taylor(erf(x), 9, 0)

3.8 曲面绘图

syms x y; ezsurf(sin(x)*sin(y), [-3*pi 3*pi -3*pi 3*pi])
ezsurf((x^2 + y^2)*cos(x^2 + y^2), [-1 1 -1 1])

3.9 动画制作

T = 0:0.01:1;
for j = 0:16
    fill(4*cos(j*pi/8) + (1/2)*cos(2*pi*T), ...
         4*sin(j*pi/8) + (1/2)*sin(2*pi*T), 'r');
    axis equal; axis([-5 5 -5 5]);
    M(j + 1) = getframe;
end
movie(M)

3.10 线性方程组求解

A1 = [3 4 5; 2 -3 7; 1 -6 1]; b = [2; -1; 3];
x = A1\b
A1*x
A2 = [3 -9 8; 2 -3 7; 1 -6 1]; b = [2; -1; 3];
x = A2\b
A3 = [1 3 -2 4; -2 3 4 -1; -4 -3 1 2; 2 3 -4 1]; b3 = [1; 1; 1; 1];
x = A3\b3
A3*x
syms a b c d x y u v;
A4 = [a b; c d]; A4\[u; v]
det(A4)

可分析矩阵的奇异性和行列式。

3.11 矩阵的秩、行列式与逆

rank(A1)
rank(A2)
rank(A3)
rank(A4)
det(A1)
inv(A1)
det(A2)
det(A3)
inv(A3)
det(A4)
inv(A4)

3.12 特征值与特征向量

[U1, R1] = eig(A1)
A1*U1 - U1*R1
[U2, R2] = eig(A2)
A2*U2 - U2*R2
[U3, R3] = eig(A3)
A3*U3 - U3*R3
[U4, R4] = eig(A4)
A4*U4 - U4*R4

可观察矩阵的特征值和特征向量的性质。

3.13 矩阵迭代

• 矩阵递推关系 ：

Xn+1 = MXn
Xn = M^nX0

• 特征值分解 ：

M = [1, 1/4, 0; 0, 1/2, 0; 0, 1/4, 1];
[U, R] = eig(M)

• 矩阵验证与极限计算 ：

M - U*R*inv(U)
Minf = U*diag([1, 1, 0])*inv(U)
Minf*X0
M^5*X0
M^10*X0

4. 实践集 C 解答

4.1 辐射强度绘图

radiation = inline(vectorize('10000/(4*pi*((x - x0)^2 + (y - y0)^2 + 1))'), 'x', 'y', 'x0', 'y0')
x = zeros(1, 5); y = zeros(1, 5);
for j = 1:5
    x(j) = 50*rand;
    y(j) = 50*rand;
end
[X, Y] = meshgrid(0:0.1:50, 0:0.1:50);
contourf(X, Y, radiation(X, Y, x(1), y(1)) + radiation(X, Y, x(2), y(2)) + radiation(X, Y, x(3), y(3)) + radiation(X, Y, x(4), y(4)) + radiation(X, Y, x(5), y(5)), 20);
colormap('gray')

从图中难以确定最佳隐藏位置，直觉上应靠近边界。

4.2 生存模拟

function r = lifeordeath(x1, y1, x0, y0)
    dosage = 10000/(4*pi*((x1 - x0)^2 + (y1 - y0)^2 + 1));
    if dosage > 50
        r = 1;
    else
        r = 0;
    end
end
x1 = 25; y1 = 25; h = 0;
for n = 1:100
    x0 = 50*rand;
    y0 = 50*rand;
    r = lifeordeath(x1, y1, x0, y0);
    h = h + r;
end
if h > 0
    disp('The Captain is dead!')
else
    disp('Picard lives!')
end

多次运行结果显示船长死亡的概率较大。

4.3 蒙特卡罗模拟

x1 = 25; y1 = 25; c = 0;
for k = 1:100
    h = 0;
    for n = 1:100
        x0 = 50*rand;
        y0 = 50*rand;
        r = lifeordeath(x1, y1, x0, y0);
        h = h + r;
    end
    if h > 0
        c = c;
    else
        c = c + 1;
    end
end
disp(['The chances of Picard surviving are = ', num2str(c/100)])

模拟得到船长的生存概率。

4.4 位置调整与生存概率

通过调整船长的位置，如 x1 = 37.5; y1 = 25 和 x1 = 50; y1 = 25 以及 x1 = 50; y1 = 50 ，再次进行蒙特卡罗模拟，可发现靠近角落时生存概率提高。

综上所述，通过调试技巧修正函数，利用各种数值计算和绘图工具解决方程求解、函数分析、矩阵运算等问题，并通过模拟实验评估生存概率，这些方法和技巧在数值计算和编程实践中具有重要的应用价值。

5. 总结与拓展应用

5.1 知识总结

• 调试与函数修正 ：掌握调试技巧对修正代码至关重要，如通过反转 x 坐标修正 shadecurves 函数，使其能正确绘制两条曲线间的区域。
• 数值计算 ：涵盖基本运算、函数计算、方程求解、因式分解、化简、幂运算、符号计算等多个方面，熟悉各种函数和命令的使用，如 solve 用于方程求解， simplify 和 simple 用于表达式化简。
• 绘图与可视化 ：可以使用 ezplot 、 plot 、 contour 、 ezsurf 等函数进行函数绘图、等高线绘图和曲面绘图，直观展示函数的性质和特点。
• 矩阵运算 ：涉及线性方程组求解、矩阵的秩、行列式、逆、特征值与特征向量的计算，以及矩阵迭代等内容，理解矩阵运算在解决实际问题中的应用。
• 模拟与概率 ：通过蒙特卡罗模拟等方法评估事件的概率，如模拟船长在辐射环境中的生存概率，了解模拟实验在实际问题中的应用。

5.2 拓展应用

5.2.1 工程领域

在工程领域，这些数值计算和编程技巧可用于电路分析、信号处理、控制系统设计等方面。例如，在电路分析中，可使用矩阵运算求解线性方程组来确定电路中的电流和电压；在信号处理中，可通过函数绘图和分析来研究信号的特性。

5.2.2 金融领域

在金融领域，可用于风险评估、投资组合优化、期权定价等问题。例如，通过蒙特卡罗模拟来评估投资组合的风险和收益，使用数值积分计算期权的价格。

5.2.3 科学研究

在科学研究中，可用于数据分析、模型求解、实验模拟等方面。例如，在物理学中，可使用导数和积分计算物理量的变化率和累积效应；在生物学中，可通过矩阵迭代模拟生物种群的演化。

5.3 操作步骤总结

以下是一些常见操作的步骤总结：
| 操作类型 | 操作步骤 |
| ---- | ---- |
| 方程求解 | 1. 使用 syms 定义符号变量；2. 使用 solve 函数输入方程和变量进行求解；3. 可通过矩阵乘法验证结果。 |
| 函数绘图 | 1. 使用 ezplot 或 plot 函数输入函数表达式和绘图范围；2. 可使用 hold on 和 hold off 进行图形叠加；3. 可使用 axis 函数调整坐标轴范围。 |
| 积分计算 | 1. 使用 int 函数输入被积函数、积分变量和积分区间进行积分计算；2. 对于不定积分，可使用 diff 函数进行验证。 |
| 矩阵运算 | 1. 定义矩阵；2. 使用 rank 计算矩阵的秩， det 计算行列式， inv 计算逆矩阵， eig 计算特征值和特征向量；3. 可使用矩阵乘法进行验证。 |
| 蒙特卡罗模拟 | 1. 定义模拟函数，如 lifeordeath ；2. 进行多次模拟实验，统计事件发生的次数；3. 计算事件发生的概率。 |

5.4 mermaid 流程图总结

以下是一些关键流程的 mermaid 流程图总结：

graph LR
    A[调试与修正代码] --> B[数值计算]
    B --> C[绘图与可视化]
    C --> D[矩阵运算]
    D --> E[模拟与概率分析]
    E --> F[拓展应用]

graph LR
    A[定义问题] --> B[选择合适的方法]
    B --> C[编写代码]
    C --> D[运行代码并调试]
    D --> E[分析结果]
    E --> F[优化和拓展]

通过以上总结和拓展，我们可以更系统地掌握数值计算和编程实践中的各种方法和技巧，并将其应用到更广泛的领域中。在实际应用中，我们可以根据具体问题选择合适的方法和工具，不断优化和拓展解决方案，以提高问题解决的效率和准确性。