《数值计算方法》学习笔记：基于MATLAB

第2章 非线性方程的求解

2.1 二分法 mbisec.m

% 程序2.1--mbisec.m
function [x, k] = mbisec(f, a, b, ep)
% 用途：用二分法求非线性方程f(x) = 0有根区间[a, b]中的一个根
% 格式：[x, k] = mbisec(f, a, b, ep) f为函数表达式，a,b为
% 区间左右端点，ep为精度，x,k分别返回近似根和二分次数
x = (a+b) / 2.0; k = 0;        
% feval:使用函数的名称或其句柄以及输入参数来计算函数的结果
while abs(feval(f, x)) > ep || (b-a > ep)
    if feval(f, x) * feval(f, a) < 0
        b = x;
    else
        a = x;
    end
    x = (a+b) / 2.0; k = k + 1;
end

% 调用格式:求方程f(x) = x - e^(-x) = 0在区间[0, 1]内的一个实根
% 取定精度ep = 10^(-5)
f = @(x)x - exp(-x);
[x, k] = mbisec(f, 0, 1, 1e-5)

2.2 迭代法 miter.m

% 程序2.2--miter.m
function [x, k] = miter(phi, x0, ep, N)
% 用途：用简单迭代法求非线性方程f(x) = 0有根区间[a, b]中的一个根
% 格式：[x, k] = miter(phi, x0, ep, N) phi为迭代函数，x0为初值，ep为精度
% (默认1e-4)，N为最大迭代次数(默认500)，x,k分别返回近似根和迭代次数
% nargin：针对当前正在执行的函数，返回函数调用中给定函数输入参数的数目。
if nargin < 4 N = 500; end
if nargin < 3 ep = 1e-4; end
k = 0;
while k < N
    x = feval(phi, x0);
    if abs(x-x0) < ep
        break;
    end
    x0 = x; k = k + 1;
end

% 调用格式:求方程f(x) = x - e^(-x) = 0在区间[0, 1]内的一个实根
% 取定精度ep = 10^(-5)，初始点为x0 = 0.5
phi = @(x)exp(-x);
[x, k] = miter(phi, 0.5, 1e-5)

% 改进:迭代函数变为"(exp(-x) + 0.5*x) / 1.5"
phi = @(x)(exp(-x) + 0.5*x) / 1.5;
[x, k] = miter(phi, 0.5, 1e-5)

2.3 Aitken加速法 maitken.m

% 程序2.3--maitken.m
function [x, k] = maitken(phi, x0, ep, N)
% 用途：用Aitken加速方法求f(x) = 0的解
% 格式：[x, k] = maitken(phi, x0, ep, N) phi为迭代函数，x0为迭代初值，
% ep为精度(默认1e-4)，N为最大迭代次数(默认500)，x,k分别返回近似根
% 和迭代次数
if nargin < 4 N = 500; end
if nargin < 3 ep = 1e-4; end
k = 0;
while k < N
    y = feval(phi, x0);
    z = feval(phi, y);
    x = x0 - (y - x0) ^ 2 / (z - 2*y + x0);
    if abs(x-x0) < ep, break; end
    x0 = x; k = k + 1;
end

% 调用格式:求方程f(x) = x - e^(-x) = 0在区间[0, 1]内的一个实根
% 取定精度ep = 10^(-5)，初始点为x0 = 0.5
phi = @(x)exp(-x);
[x, k] = maitken(phi, 0.5, 1e-5)

2.4 牛顿法 mnewton.m

% 程序2.4--mnewton.m
function [x, k] = mnewton(f, df, x0, ep, N)
% 用途：用牛顿法求解非线性方程f(x) = 0
% 格式：[x, k] = maitken(phi, x0, ep, N) f和df分别为表示f(x)
% 及其导数，x0为迭代初值，ep为精度(默认1e-4)，N为最大迭
% 代次数(默认500)，x,k分别返回近似根和迭代次数
if nargin < 5, N = 500; end
if nargin < 4,ep = 1e-4; end
k = 0;
while k < N
    x = x0 - feval(f, x0) / feval(df, x0);
    if abs(x-x0) < ep
        break;
    end
    x0 = x; k = k + 1;
end

% 调用格式:求方程f(x) = x - e^(-x) = 0在区间[0, 1]内的一个实根
% 取定精度ep = 10^(-5)，初始点为x0 = 0.5
f = @(x)x - exp(-x);
df = @(x)1 + exp(-x);
[x, k] = mnewton(f, df, 0.5, 1e-5)

2.5 割线法 mqnewt.m

% 程序2.5--mqnewt.m
function [x, k] = mqnewt(f, x0, x1, ep, N)
% 用途：用割线法求解非线性方程f(x) = 0
% 格式：[x, k] = mqnewt(f, x0, x1, ep, N) f为f(x)的表达式，
% x0,x1为迭代初值，ep为精度(默认1e-4)，N为最大迭代次
% 数(默认500)，x,k分别返回近似根和迭代次数
if nargin < 5, N = 500; end
if nargin < 4, ep = 1e-4; end
k = 0;
while k < N
    x = x1 - (x1-x0) * feval(f, x1) / (feval(f, x1) - feval(f, x0));
    if abs(x-x1) < ep
        break;
    end
    x0 = x1; x1 = x;
    k = k + 1;
end

% 调用格式:求方程f(x) = e^x + x^2 - 2 = 0的实根
% 取定精度ep = 10^(-5)，初始点为x0 = 0.0,x1 = 1.0
f = @(x)exp(x) + x^2 - 2;
[x, k] = mqnewt(f, 0.0, 1.0, 1e-5)

第3章 线性方程组的直接解法

3.1 顺序高斯消去法 mgauss.m

% 程序3.1--mgauss.m
function [x] = mgauss(A, b, flag)
% 用途：顺序高斯消去法解线性方程组Ax = b
% 格式：[x] = maguss(A, b, flag), A为系数矩阵，b为右端项，若flag = 0,
% 则不显示中间过程，否则显示中间过程，默认为0，x为解向量
if nargin < 3, flag = 0; end
n = length(b);
% 消元过程
for k = 1 : (n - 1)
    m = A(k+1 : n, k) / A(k, k);
    A(k+1 : n, k+1 : n) = A(k+1 : n, k+1 : n) - m * A(k, k+1 : n);
    b(k+1 : n) = b(k+1 : n) - m * b(k);
    A(k+1 : n, k) = zeros(n-k, 1);
    if flag ~= 0, Ab = [A, b], end
end
% 回代过程
x = zeros(n, 1);
x(n) = b(n) / A(n, n);
for k = n-1 : -1 : 1
    x(k) = (b(k) - A(k, k+1 : n) * x(k+1 : n)) / A(k, k);
end

%调用格式
A = [1 1 1 1; -1 2 -3 1; 3 -3 6 -2; -4 5 2 -3];
b=[10 -2 7 0]';
x = mgauss(A, b)

3.2 列主元高斯消去法 mgauss2.m

% 程序3.2--mgauss2.m
function [x] = mgauss2(A, b, flag)
% 用途：列主元高斯消去法解线性方程组Ax = b
% 格式：[x] = maguss(A, b, flag),A为系数矩阵，b为右端项，若
% flag = 0(默认),则不显示中间过程，否则显示中间过程，x为解向量
if nargin < 3, flag = 0; end
n = length(b);
for k = 1 : (n-1) %选主元
    [ap, p] = max(abs(A(k : n, k)));
    p = p + k - 1;
    if p > k
        A([k p], :) = A([p k], :);
        b([k p], :) = b([p k], :);
    end
    % 消元
    m = A(k+1 : n, k) / A(k, k);
    A(k+1 : n, k+1 : n) = A(k+1 : n, k+1 : n) - m * A(k, k+1 : n);
    b(k+1 : n) = b(k+1 : n) - m * b(k);
    A(k+1 : n, k) = zeros(n-k, 1);
    if flag ~= 0,Ab = [A, b], end
end
% 回代
x = zeros(n, 1);
x(n) = b(n) / A(n, n);
for k = n-1 : -1 : 1
    x(k) = (b(k) - A(k, k+1 : n) * x(k+1 : n)) / A(k, k);
end

%调用格式
A=[2 -1 4 -3 1; -1 1 2 1 3; 4 2 3 3 -1; -3 1 3 2 4; 1 3 -1 4 4];
b=[11 14 4 16 18]';
x = mgauss2(A, b); x = x'

3.3 LU分解法 mlu.m

% 程序3.3--mlu.m
function [x, L, U] = mlu(A, b)
% 用途：用LU分解法解方程组Ax = b
% 格式：[x, L, U] = mlu(A, b) A为系数矩阵，b为右端向量
% x返回解向量，L返回下三角矩阵，U返回上三角矩阵
n = length(b);
% LU分解
U = zeros(n, n); L = eye(n, n);
U(1, :) = A(1, :); L(2:n, 1) = A(2:n, 1) / U(1, 1);
for k = 2 : n
    U(k, k:n) = A(k, k:n) - L(k, 1:k-1) * U(1:k-1, k:n);
    L(k+1:n, k) = (A(k+1:n, k)- L(k+1:n, 1:k-1) * U(1:k-1, k)) / U(k, k);
end
% 解下三角方程组Ly = b
y = zeros(n, 1); y(1) = b(1);
for k = 2 : n
    y(k) = b(k) - L(k, 1:k-1) * y(1 : k-1);
end
% 解上三角方程组 Ux = y
x = zeros(n, 1); x(n) = y(n) / U(n, n);
for k = n-1 : -1 : 1
    x(k) = (y(k) - U(k, k+1:n) * x(k+1 : n)) / U(k, k);
end

% 调用格式
A = [2 -1 4 -3 1; -1 1 2 1 3; 4 2 3 3 -1; -3 1 3 2 4; 1 3 -1 4 4];
b = [11 14 4 16 18]';
[x, L, U] = mlu(A, b)

3.4 列主元LU分解法 mzlu.m

% 程序3.4--mzlu.m
function [x, L, U, P] = mzlu(A, b)
% 用途：用列主元LU分解法解方程组Ax = b
% 格式：[x, L, U, P] = mzlu(A, b) A为系数矩阵，b为右端向量，
% x返回解向量，L返回下单位三角矩阵，U返回上三角矩阵
% P返回选主元时记录行交换的置换阵
n = length(b);
P = eye(n); % P记录选择主元时候所进行的行变换
% 列主元LU分解
for k = 1 : n
    A(k:n, k) = A(k:n, k) - A(k:n, 1:k-1) * A(1:k-1, k);
    [s, m] = max(abs(A(k:n, k))); % 选列主元
    m = m + k - 1;
    if m ~= k
       A([k m], :) = A([m k], :);
       P([k m], :) = P([m k], :);
       % b([k m],:) = b([m k], :);
    end
    A(k+1:n, k) = A(k+1:n, k) / A(k, k);
    A(k, k+1:n) = A(k, k+1:n) - A(k, 1:k-1) * A(1:k-1, k+1:n);
end
L = tril(A, -1) + eye(n, n);
U = triu(A);
% 解下三角矩阵Ly = b
newb = P * b;
% newb = b;
y = zeros(n, 1);
for k = 1 : n
    j = 1 : k - 1;
    y(k) = (newb(k) - L(k,j) * y(j)) / L(k, k);
end
% 解上三角方程组Ux = y
x = zeros(n, 1);
for k = n : -1 : 1
    j = k+1 : n;
    x(k) = (y(k) - U(k, j) * x(j)) / U(k, k);
end

%调用格式
A = [2 -1 4 -3 1; -1 1 2 1 3; 4 2 3 3 -1; -3 1 3 2 4; 1 3 -1 4 4];
b = [11 14 4 16 18]';
[x, L, U, P] = mzlu(A, b)

3.5 乔列斯基分解法 mchol.m

% 程序3.5--mchol.m
function [x, L, D] = mchol(A, b)
% 用途：用乔列斯基分解法解对称正定方程组Ax = b
% LDL'分解
n = length(b); D = zeros(1, n); L = eye(n, n);
U(1, :) = A(1, :);
L(2:n, 1) = A(2:n, 1) / U(1, 1);
for k = 2 : n
    U(k, k:n) = A(k, k:n) - L(k, 1:k-1) * U(1:k-1, k:n);
    L(k+1:n, k) = U(k, k+1:n) / U(k, k);
end
D = diag(diag(U));
% 求解下三角方程组Ly = b(向前消去法)
y = zeros(n, 1);
y(1) = b(1);
for k = 2 : n
    y(k) = b(k) - L(k, 1:k-1) * y([1:k-1]);
end
% 求解对角方程组Dz = y
for k = 1 : n
    z(k) = y(k) / D(k, k);
end
% 求解上三角方程组L'x = z(回代法)
x = zeros(n, 1);
U = L';
x(n) = z(n);
for k = n-1 : -1 : 1
    x(k) = z(k) -U(k, k+1:n) * x(k+1:n);
end

% 调用格式
A = [4 2 -2; 2 2 -3; -2 -3 14];
b = [4 1 0]';
[x, L, D] = mchol(A, b)

3.6 追赶法 mchase.m

% 程序3.6--mchase.m
function [x] = mchase(a, b, c, d)
% 用途：追赶法解三对角方程组Ax = d
% 格式：[x] = mchase(a, b, c, d),a为次下对角线元素向量，b主对角
% 元素向量，c为次上对角线元素向量，d为右端向量，x返回解向量
n = length(b);
for k = 2 : n
    b(k) = b(k) - a(k) / b(k-1) * c(k-1);
    d(k) = d(k) - a(k) / b(k-1) * d(k-1);
end
x(n) = d(n) / b(n);
for k = n-1 : -1 : 1
    x(k) = (d(k) - c(k) * x(k+1)) / b(k);
end

% 调用格式
b = [2 3 4 5]';
a1 = 0; c4 = 0;
a = [a1 -1 -2 -3]';
c = [-1 -2 -3 c4]';
d = [6 1 -2 1]';
[x] = mchase(a, b, c, d)

第4章 线性方程组的迭代解法

4.1 雅可比迭代法 mjacobo.m

% 程序4.1--mjacobo.m
function [x, iter] = mjacobi(A, b, x, ep, N)
% 用途：用雅可比迭代法解线性方程组Ax = b
% 格式：[x, iter] = mjacobi(A, b, x, ep, N) A为系数矩阵，b为右端向
% 量，x为初始向量(默认零向量)，ep为精度(默认1e-6)，N为最大迭
% 代次数(默认500次)，返回参数x，iter分别为近似解向量和迭代次数
if nargin < 5, N = 500; end
if nargin < 4, ep = 1e-6; end
if nargin < 3, x = zeros(size(b)); end
D = diag(diag(A));
for iter = 1 : N
    x = D \ ((D-A) * x + b);
    err = norm(b - A*x) / norm(b);
    if err < ep, break; end
end

% 测试用例
A = [0.76 -0.01 -0.14 -0.16;
    -0.01 0.88 -0.03 0.05;
    -0.14 -0.03 1.01 -0.12;
    -0.16 0.05 -0.12 0.72];
b = [0.68 1.18 0.12 0.74]';
[x, iter] = mjacobi(A, b)

4.2 高斯-塞德尔迭代法 mseidel.m

% 程序4.2--mseidel.m
function [x, iter] = mseidel(A, b, x, ep, N)
% 用途：用高斯-塞德尔迭代法解线性方程组Ax = b
% 格式：[x, iter] = mseidel(A, b, x, ep, N) A为系数矩阵，b为右端向
% 量，x为初始向量(默认零向量)，ep为精度(默认1e-6)，N为最大迭
% 代次数(默认500次)，返回参数x，iter分别为近似解向量和迭代次数
if nargin < 5, N = 500; end
if nargin < 4,ep = 1e-6; end
if nargin < 3, x = zeros(size(b)); end
D = diag(diag(A)); L = D - tril(A); U = D - triu(A);
for iter = 1 : N
    x = (D - L) \ (U * x + b);
    err = norm(b - A*x) / norm(b);
    if err < ep, break; end
end

% 测试用例
A = [0.76 -0.01 -0.14 -0.16;
    -0.01 0.88 -0.03 0.05;
    -0.14 -0.03 1.01 -0.12;
    -0.16 0.05 -0.12 0.72];
b = [0.68 1.18 0.12 0.74]';
[x, iter] = mseidel(A, b)

4.3 SOR迭代法 msor.m

% 程序4.3--msor.m
function [x, iter] = msor(A, b, omega, x, ep, N)
% 用途：用SOR迭代法解线性方程组Ax = b
% 格式：[x, iter] = msor(A, b, omega, x, ep, N)其中A为系数矩阵，b为右端向量，
% omega为松弛因子(默认1.2)，x为初始向量(默认零向量)，ep为精度(默认1e-6)，
% N为最大迭代次数(默认500次)，返回参数x，iter分别为近似解向量和迭代次数
if nargin < 6, N = 500; end
if nargin < 5, ep = 1e-6; end
if nargin < 4,x = zeros(size(b)); end
if nargin < 3, omega = 1.2; end
D = diag(diag(A)); L = D - tril(A); U = D - triu(A);
for iter = 1 : N
    x = (D - omega * L) \ (((1-omega) * D + omega * U) * x + omega * b);
    err = norm(b - A*x) / norm(b);
    if err < ep, break; end
end

% 测试用例
A = [0.76 -0.01 -0.14 -0.16;
    -0.01 0.88 -0.03 0.05;
    -0.14 -0.03 1.01 -0.12;
    -0.16 0.05 -0.12 0.72];
b = [0.68 1.18 0.12 0.74]';
[x, iter] = msor(A, b, 1.05)

4.4 最速下降法 mgrad.m

% 程序4.4--mgrad.m
function [x, iter] = mgrad(A, b, x, ep, N)
% 用途：用最速下降法解线性方程组Ax = b
% 格式：[x, iter] = mgrad(A, b, x, ep, N) 其中A为系数矩阵，b为右端向
% 量，x为初始向量(默认零向量)，ep为精度(默认1e-6)，N为最大迭代次
% 数(默认500次)，返回参数x，iter分别为近似解向量和迭代次数
if nargin < 5, N = 1000; end
if nargin < 4, ep = 1e-6; end
if nargin < 3, x = zeros(size(b)); end
for iter = 1 : N
    r = b - A * x;
    if norm(r) < ep, break; end
    a = r' * r / (r' * A * r);
    x = x + a * r;
end

% 测试用例
A = [0.76 -0.01 -0.14 -0.16;
    -0.01 0.88 -0.03 0.05;
    -0.14 -0.03 1.01 -0.12;
    -0.16 0.05 -0.12 0.72];
b = [0.68 1.18 0.12 0.74]';
[x, iter] = mgrad(A,b)

4.5 共轭梯度法 mcg.m

% 程序4.5--mcg.m
function [x, iter] = mcg(A, b, x, ep, N)
% 用途：用共轭梯度法解线性方程组Ax = b
% 格式：[x, iter] = mcg(A, b, x, ep, N) 其中A为系数矩阵，b为右端
% 向量，x为初始向量(默认零向量)，ep为精度(默认1e-6)，N为最大迭
% 代次数(默认500次)，返回参数x，iter分别为近似解向量和迭代次数
if nargin < 5, N = 1000; end
if nargin < 4, ep = 1e-6; end
if nargin < 3, x = zeros(size(b)); end
r = b - A * x;
for iter = 1 : N
    rr = r' * r;
    if iter == 1
        p = r;
    else
        beta = rr / rr1;
        p = r + beta * p;
    end
    q = A * p;
    alpha = rr / (p' * q);
    x = x+ alpha * p;
    r = r- alpha * q;
    rr1 = rr;
    if norm(r) < ep, break; end
end

% 测试用例
A = [0.76 -0.01 -0.14 -0.16;
    -0.01 0.88 -0.03 0.05;
    -0.14 -0.03 1.01 -0.12;
    -0.16 0.05 -0.12 0.72];
b = [0.68 1.18 0.12 0.74]';
[x, iter] = mcg(A, b)

4.6 广义极小残量法 mgmres.m

% 程序4.6--mgmres.m
function [x, iter] = mgmres(A, b, x, ep, N)
% 用途：用广义极小残量法解线性方程组Ax = b
% 格式：[x, iter] = mgmres(A, b, x, ep, N) 其中A为系数矩阵，b为右端
% 向量，x为初始向量(默认零向量)，ep为精度(默认1e-6)，N为最大迭
% 代次数(默认500次)，返回参数x，iter分别为近似解向量和迭代次数
if nargin < 5, N = 500; end
if nargin < 4, ep = 1e-6; end
if nargin < 3, x = zeros(size(b)); end
n = length(b);
V(1:n, 1:n+1) = zeros(n, n+1);
H(1:n+1, 1:n) = zeros(n+1, n);
cs(1:n) = zeros(n, 1);
sn(1:n) = zeros(n, 1);
e1 = zeros(n, 1); e1(1) = 1.0;
for iter = 1 : N
    r = b - A * x;
    if norm(r) < ep, break; end
    beta = norm(b); V(:,1) = r / beta;
    s = beta * e1;
    for i = 1 : n
        u = A * V(:, i);
        for k = 1 : i
            H(k, i) = u' * V(:, k);
            u = u - H(k, i) * V(:, k);
        end
        H(i+1, i) = norm(u);
        V(:, i+1) = u / H(i+1, i);
        for k = 1 : i-1
            temp = cs(k) * H(k, i) + sn(k) * H(k+1, i);
            H(k+1, i) = -sn(k) * H(k, i) + cs(k) * H(k+1, i);
            H(k, i) = temp;
        end
        cs(i) = H(i, i) / sqrt(H(i, i) ^ 2 + H(i+1, i) ^ 2);
        sn(i) = H(i+1, i) / sqrt(H(i, i) ^ 2 + H(i+1, i) ^ 2);
        temp = cs(i) * s(i); s(i+1) = -sn(i) * s(i); s(i) = temp;
        H(i, i) = cs(i) * H(i, i) + sn(i) * H(i+1, i);
        H(i+1, i) = 0;
        if abs(s(i+1)) < ep
            y = H(1:i, 1:i) \ s(1:i);
            x = x + V(:, 1:i) * y;
            break;
        end
    end
end

% 测试用例
A = [0.76 -0.01 -0.14 -0.16;
    -0.01 0.88 -0.03 0.05;
    -0.14 -0.03 1.01 -0.12;
    -0.16 0.05 -0.12 0.72];
b = [0.68 1.18 0.12 0.74]';
[x, iter] = mgmres(A, b)        

第5章 插值法和最小二乘拟合

5.1 拉格朗日插值法 mlagr.m

% 程序5.1--mlagr.m
function yp = mlagr(x, y, xp)
% 用途：拉格朗日插值法求解
% 格式：yp = mlagr(x, y, xp)，x是节点向量，y是节点对应的函
% 数值向量，xp是插值点(可以是多个)，yp返回插值结果
n = length(x); m = length(xp);
yp = zeros(1, m); c1 = ones(n-1, 1); c2 = ones(1, m);
for i = 1 : n
    xb = x([1:i-1, i+1:n]);
    yp = yp + y(i) * prod((c1 * xp - xb' * c2) ./ (x(i) - xb' * c2));
end

% 测试用例
x = pi * [1/6 1/4];
y = [0.5 0.7071]; xx = 2 * pi / 9;
yy1 = mlagr(x, y, xx)

5.2 牛顿插值法 mnewp.m

% 程序5.2--mnewp.m(存在问题)
function yy = mnewp(x, y, xx)
% 用途：牛顿插值法求解
% 格式：yy = mnewp(x, y, xx)，x是节点向量，y是节点对应的函
% 数值向量，xx是插值点(可以是多个)，yy返回插值结果
n = length(x);
syms t; yy = y(1);
y1 = 0; lx = 1;
for i = 1 : n-1
    for j = i+1 : n
        y1(j) = (y(j-1) - y(j)) / (x(j-i) - x(j)); %计算差商
    end
    c(i) = y1(i+1); lx = lx * (t - x(i));
    yy = yy + c(i) * lx;   %计算牛顿插值多项式的值
    y = y1;
end
if nargin == 3
    yy = subs(yy, 't', xx);
else
    yy = collect(yy);
    yy = vpa(yy,6);
end

% 测试用例
x = pi * [1/6 1/4 1/3 1/2];
y = [0.5 0.7071 0.8660 1];
xx = 2 * pi / 9;
yy = mnewp(x, y, xx)

5.3 Runge现象 mrunge.m

% 程序5.3--mrunge.m(存在问题)
function mrunge()
% 高阶插值的Runge现象
xx = -5 : 0.05 : 5; y = 1 ./ (1 + xx .^ 2);
x1 = -5 : 2 : 5; y1 = 1 ./ (1 + x1 .^ 2);
x2 = -5 : 1 : 5; y2 = 1./ (1 + x2 .^ 2);
yy1 = mlagr(x1, y1, xx);
yy2 = mlagr(x2, y2, xx);
plot(xx, yy1, 'k-.'); hold on
plot(xx, yy2, 'k-.'); plot(xx, y, 'k');
legend('10次多项式插值', '5次多项式插值', '被插函数的图形');
axis([-5 5 -0.5 2])

5.4 分段线性插值 mpiecel.m

% 程序5.4--mpiecel.m
function yy = mpiecel(x, y, xx)
% 用途：分段线性插值
% 格式：yy = mpiecel(x, y, xx)，x是节点向量，y是节点对应的函
% 数值向量，xx是插值点(可以是多个)，yy返回插值结果
n = length(x);
for j = 1 : length(xx)
	for i = 2 : n
		if xx(j) < x(i)
			yy(j) = y(i-1) * (xx(j)-x(i)) / (x(i-1)-x(i)) + y(i) * (xx(j)-x(i-1)) / (x(i)-x(i-1));
			break;
		end
	end
end

% 测试用例
x = pi * [1/6 1/4 1/3 1/2];
y = [0.5 0.7071 0.8660 1];
xx = 2 * pi / 9;
yy = mnewp(x, y, xx)

5.5 三次样条插值 mspline.m

% 程序5.5--mspline.m
function m = mspline(x, y, dy0, dyn, xx)
% 用途：三次样条插值(一阶导数边界条件)
% 格式：m = mspline(x, y, dy0, dyn, xx)，
% x，y分别为n个节点的横坐标所组成的向量及纵坐标所组成的向量，dy0,dyn在左右两端
% 点的一阶导数，如果xx缺省，则输出各节点的一阶导数值，否则，m为xx的三次样条插值
n = length(x) - 1; % 计算小区间的个数
h = diff(x); lambda = h(2:n) ./ (h(1:n-1) + h(2:n)); mu = 1 - lambda;
theta = 3 * (lambda .* diff(y(1:n)) ./ h(1:n-1) + mu .* diff(y(2:n+1)) ./ h(2:n));
theta(1) = theta(1) - lambda(1) * dy0;
theta(n-1) = theta(n-1) - lambda(n-1) * dyn;
% 追赶法解三对角方程组
dy = machase(lambda, 2*ones(1:n-1), mu, theta);
% 若给出插值点，计算相应的插值
m = [dy0; dy; dyn];
if nargin >= 5
    s = zeros(size(xx));
    for i = 1 : n
        if i == 1
            kk = find(xx <= x(2));
        elseif i == n
            kk = find(xx > x(n));
        else
            kk = find(xx > x(i) & xx <= x(i+1));
        end
        xbar = (xx(kk) - x(i)) / h(i);
        s(kk) = alpha0(xbar) * y(i) + alpha1(xbar) * y(i+1)+...
            +h(i) * beta0(xbar) * m(i) + h(i) * beta1(xbar) * m(i+1);
    end
    m = s;
end
% 追赶法
function x = machase(a, b, c, d)
n = length(a);
for k = 2 : n
    b(k) = b(k) - a(k) / b(k-1) * c(k-1);
    d(k) = d(k)- a(k) / b(k-1) * d(k-1);
end
x(n) = d(n) / b(n);
for k = n-1 : -1 : 1
    x(k) = (d(k) - c(k) * x(k+1)) / b(k);
end
x = x(:);
% 基函数
function y = alpha0(x)
y = 2 * x .^ 3 - 3 * x .^ 2 + 1;
function y = alpha1(x)
y = -2 * x .^ 3 + 3 * x .^ 2;
function y = beta0(x)
y = x .^ 3 - 2 * x .^ 2 + x;
function y = beta1(x)
y = x .^ 3 -x .^ 2;

% 测试用例
x = [-1 0 1 2]; y = [-1 0 1 0];
xx = [-0.8 -0.3 0.2 0.7 1.2 1.7];
yy = mspline(x, y, 0, -1, xx)

5.6 多项式拟合 mpfit.m

% 程序5.6--mpfit.m
function p = mpfit(x, y, m)
% 用途：多项式拟合
% 格式：p = mpfit(x, y, m)
% x,y为数据向量，m为拟合多项式次数，p返回多项式系数降幂排列
A = zeros(m+1, m+1);
for i = 0 : m
    for j = 0 : m
        A(i+1, j+1) = sum(x .^ (i + j));
    end
    b(i+1) = sum(x .^ i .* y);
end
a = A \ b';
p = fliplr(a');%按降幂排列

% 测试用例
x = -2:2; y = [-1 -1 0 1 1];
p = mpfit(x, y, 3)

第6章 数值积分和数值微分

6.1 复化中点公式 mintm.m

% 程序6.1--mintm.m
function s = mintm(f, a, b, n)
% 用途:用复化中点公式求积分.
% 格式:s = mintm(f, a, b, n) f为被积函数,a,b为积分
% 区间的左右端点,n为区间的等分数,s返回积分近似值
h = (b - a) / n;
x = linspace(a + h / 2, b - h / 2, n);
y = feval(f, x);
s = h * sum(y)

% 测试用例
format long
f = @(x)4 ./ (1 + x .^ 2);
s = mintm(f, 0, 1, 20)

6.2 复化梯形公式 mtrap.m

% 程序6.2--mtrap.m
function s = mtrap(f, a, b, n)
% 用途:用复化梯形公式求积分.
% 格式:s = mtrap(f, a, b, n) f为被积函数,a,b为积分
% 区间的左右端点,n为区间的等分数,s返回积分近似值
h = (b - a) / n;
x = linspace(a, b, n+1);
y = feval(f, x);
s = 0.5 * h * (y(1) + 2 * sum(y(2:n)) + y(n+1));

% 测试用例
format long
f = @(x)4 ./ (1 + x .^ 2);
s = mtrap(f, 0, 1, 20)

6.3 复化辛普森公式 msimp.m

% 程序6.3--msimp.m
function s = msimp(f, a, b, n)
% 用途:用复化辛普森公式求积分.
% 格式: s = msimp(f, a, b, n) f为被积函数;a,b为积分
% 区间的左右端点,n为区间的等分数,s返回积分近似值
h = (b - a) / n;
x = linspace(a, b, 2 * n + 1);
y = feval(f, x);
s = (h/6) * (y(1) + 2 * sum(y(3:2:2*n-1)) + 4 * sum(y(2:2:2*n)) + y(n+1));

% 测试用例
format long
f = @(x)4 ./ (1 + x .^ 2);
s = msimp(f, 0, 1, 20)

6.4 复化梯形公式 mvtrap.m

% 程序6.4--mvtrap.m
function [T2, k, n] = mvtrap(f, a, b, eps)
% 用途:用复化梯形公式求积分.
% 格式: [T2, k] = mvtrap(f, a, b, eps) f为被积函数,a,b为积分区间的左右端
% 点,eps为控制精度,T2返回积分近似值，k为对分区间的次数,n为区间的等分数
h = b - a;
T1 = h * (feval(f, a) + feval(f, b)) / 2;
T2 = T1 / 2 + (h/2) * feval(f, a + h / 2);
k = 1; n = 2;
while(abs(T2-T1) > eps)
    h = h / 2; T1 = T2;
    T2 = T1 / 2 + (h/2) * sum(feval(f, a + h / 2 : h : b - h / 2));
    k = k + 1; n = 2 * n;
end

% 测试用例
format long
f = @(x)4 ./ (1 + x .^ 2);
[T, k, n] = mvtrap(f, 0, 1, 1.e-10)

6.5 龙贝格公式 mromb.m

% 程序6.5--mromb.m
function [T, n] = mromb(f, a, b, eps)
% 用途:用龙贝格公式求积分
% 格式: [R, n] = mromb(f, a, b, eps), f是被积函數，[a, b]是积分
% 区间,eps控制精度，R返回积分近似值,n返回区间等分数
if nargin < 4, eps = 1e-6; end
h = b - a;
R(1,1) = (h/2) * (feval(f, a) + feval(f, b));
n = 1; J = 0; err = 1;
while (err > eps)
    J = J + 1; h = h / 2; S = 0;
    for i = 1 : n
        x = a + h * (2 * i - 1);
        S = S + feval(f, x);
    end
    R(J+1, 1) = R(J, 1) / 2 + h * S;
    for k = 1 : J
        R(J+1, k+1) = (4^k * R(J+1, k) - R(J, k)) / (4^k - 1);
    end
    err = abs(R(J+1, J+1) - R(J+1, J));
    n = 2 * n;
end
R; %龙贝格表
T = R(J+1, J+1);

% 测试用例
format long
f = @(x)4 ./ (1 + x.^2);
[T, n] = mromb(f, 0, 1, 1.e-10)

6.6 定步长高斯求积公式 mgsint.m

% 程序6.6--mgsint.m
function g = mgsint(f, a, b, n, m)
% 用途:用定步长高斯求积公式求函数的积分
% 格式: g = mgsint(f, a, b, n, m) fname是被积函数,
% a,b分别为积分下上限,n为等分数,m为每段高斯点数
switch m
    case 1
        t = 0; A = 1;
    case 2
        t = [-1 / sqrt(3), 1 / sqrt(3)]; A = [1, 1];
    case 3
        t = [-sqrt(0.6), 0.0, sqrt(0.6)]; A = [5/9, 8/9, 5/9];
    case 4
        t = [-0.861136, -0.339981, 0.339981, 0.861136];
        A = [0.347855, 0.652145, 0.652145, 0.347855];
    case 5
        t = [-0.906180, -0.538469, 0.0, 0. 538469, 0.906180];
        A = [0.236927, 0.478629, 0.568889, 0.478629, 0.236927];
    case 6
        t = [-0.932470, -0.661209, -0.238619, 0.238619, 0.661209, 0.932470];
        A = [0.171325, 0.360762, 0.467914, 0.467914, 0.360762, 0.171325];
    otherwise
        error('本程序高斯点数只能取1,2,3,4,5,6!');
end
x = linspace(a, b, n+1); g = 0;
for i = 1 : n
    g = g + gsint(f, x(i), x(i+1), A, t);
end
% 子函数
function g = gsint(f, a, b, A, t)
g = (b-a) / 2 * sum(A .* feval(f, (b-a)/2*t + (a+b)/2));

% 测试用例
format long
f = @(x)4 ./ (1 + x .^ 2);
g = mgsint(f, 0, 1, 2, 3)
g = mgsint(f, 0, 1, 4, 4)
g = mgsint(f, eps, 1, 2, 3)
g = mgsint(f, eps, 1,,3)

第7章 矩阵特征值问题的数值方法

7.1 乘幂法 meigpower.m

% 程序7.1--meigpower.m
function [lam, v, k] = meigpower(A, x, eps, N)
% 用途:用乘幂法求矩阵的模最大特征值和对应的特征向量
% 格式: [lam, v, k] = meigpower(A, xO, eps, N)
% A为n阶方阵，x为初始向量，eps控制精度，N为最大迭代次数.
% 1am返回按模最大的特征值，v返回对应的特征向量，k返回迭代次数.
if nargin < 4, N = 500; end
if nargin < 3, eps = 1e-6; end
m = 0; k = 0; err = 1;
while(err > eps)
    v = A * x;
    [m1, t] = max(abs(v));
    m1 = v(t);
    x = v / m1;
    err = abs(m1 - m);
    m = m1;
    k = k + 1;
end
lam = m1;
v = x;

% 测试用例
A = [-1 2 3; 2 -3 5; 3 5 -2];
x = [1 1 1]';
[lam, v, k] = meigpower(A, x)

7.2 原点位移乘幂法 meigpowerp.m

% 程序7.2--meigpowerp.m
function [lam, v, k] = meigpowerp(A, x, alpha, eps ,N)
% 用途:用原点位移乘幂法求矩阵的模最大特征值和对应的特征向量
% 格式: [1am, v, k] = mapowerp(A, x, alpha, eps, N)
% A为n阶方阵，x为初始向量，eps为控制精度，N最大迭代次数，alpba为原点位
% 移参数.1am返回按模最大的特征值，v返回对应的特征向量，k返回迭代次数.
if nargin < 5, N = 500; end
if nargin < 4, eps = 1e-6; end
m = 0; k = 0; err = 1;
A = A - alpha * eye(length(x));
while((k < N) & (err > eps))
    v = A * x; [m1,t] = max(abs(v));
    m1 = v(t); x = v / m1;
    err = abs(m1 - m) ;
    m = m1; k = k + 1;
end
lam = m1 + alpha; v = x;

% 测试用例
A = [-1 2 3; 2 -3 5; 3 5 -2];
x = [1 1 1]';
[lam, v, k] = meigpowerp(A, x, 2)

7.3 原点位移加速反幂法 meiginvpower.m

% 程序7.3--meiginvpower.m
function [lam, v, k] = meiginvpower(A, x, alpha, eps, N)
% 用途:用原点位移加速反幂法求矩阵的模最大特征值和对应的特征向量
% 格式: [1am, V, k] = meiginvpower(A, x, alpha, eps, N)
% A为n阶方阵,x为初始向量,eps为控制精度,N为最大迭代次数,alpha为模最大的
% 近似特征值.lam返回按模最大的特征值,v返回对应的特征向量,k返回迭代次数
if nargin < 5, N = 500; end
if nargin < 4, eps = 1e-6; end
m = 0.5; k = 0; err = 1;
A = A - alpha * eye(length(x));
[L, U, P] = lu(A);
while(k < N) & (err > eps)
    [m1, t] = max(abs(x));
    m1 = x(t); v = x / m1;
    z = L \ (P * v); x = U \ z;
    err = abs(1 / m1 - 1 / m);
    k = k + 1; m = m1;
end
lam = alpha + 1 / m;

% 测试用例
A = [-1 2 3; 2 -3 5; 3 5 -2];
x = [1 1 1]';
[lam, v1, k1] = meiginvpower(A, x, -7.5)
[lam, v2, k2] = meiginvpower(A, x, 4.5)
[lam, v3, k3] = meiginvpower(A, x, -3.0)

7.4 Jacobi迭代法 meigjacobi.m

% 程序7.4--meigjacobi.m
function [D, V] = meigjacobi(A, eps)
% 用途:用Jacobi迭代法求实对称矩阵A的特征值和特征向量
% 格式: [V,D] = meigjacobi(A, eps),
% A为n阶对称方阵，eps为容许误差，V返回特征向量矩阵，
% D是n阶对角矩阵，其对角元为矩阵A的n个特征值
if nargin < 2, eps = 1e-6; end
[n, n] = size(A); D = [ ]; V = eye(n);
% 计算A的非对角元绝对值最大元素所在的行p和列q
[w1, p] = max(abs(A - diag(diag(A))));
[w2, q] = max(w1); p = p(q);
while(1)
    d = (A(q,q) - A(p,p)) / (2 * A(p,q));
    if(d > 0)
        t = -d + sqrt(d^2 + 1);
    else if(d < 0)
        t = -d - sqrt(d^2 + 1);
        else
            t = 1;
        end
    end
    c = 1 / sqrt(t^2 + 1); s = c * t;
    R = [c s;-s c];
    A([p q], :) = R' * A([p q], :);
    A(:, [p q]) = A(:, [p q]) * R;
    V(:, [p q]) = V(:, [p q]) * R;
    [w1, p] = max(abs(A - diag(diag(A))));
    [w2, q] = max(w1); p = p(q);
    if (abs(A(p, q)) < eps * sqrt(sum(diag(A) .^ 2) / n))
        break;
    end
end
D = diag(diag(A));

% 测试用例
A = [-1 2 3; 2 -3 5; 3 5 -2];
[D, V] = meigjacobi(A)

7.5 构造Householder变换矩阵H mhouseh.m

% 程序7.5--mhouseh.m
function H = mhouseh(x)
% 用途:对于向量x，构造Householder变换矩阵H,使得Hx=(*.,....,0)'
% 格式: function H = mhouseh(x)
% x为输入列向量，H返回Householder变换矩阵
n = length(x); I = eye(n); sn = sign(x(1));
if sn == 0, sn = 1; end
z = x(2:n);
if (norm(z, inf) == 0), H = I; return; end
a = sn * norm(x, 2);
u = x; u(1) = u(1) + a;
rho = a * (a + x(1));
H = I - 1.0 / rho * u * u';

% 测试用例
x = [2 1 -3 4]';
H = mhouseh(x); x1 = H * x

7.6 Householder变换化矩阵A为上Hessenberg矩阵 mhessen.m

% 程序7.6--mhessen.m
function A = mhessen(A)
% 用途:用Householder变换化矩阵A为上Hessenberg矩阵.
% 输入: n阶实方阵A
% 输出: A的上Hessenberg形
% 调用函数: mhouseh.m
[n, n] = size(A);
for k = 1 : (n-2)
    x = A(k+1:n, k); H = mhouseh(x);
    A(k+1:n, 1:n) = H * A(k+1:n, 1:n);
    A(1:n, k+1:n) = A(1:n, k+1:n) * H;
end

% 测试用例
A = [-1 2 3 5; 2 -3 8 1; 3 8 -2 7; 5 1 7 6];
A = mhessen(A)

7.7 Givens变换对上Hessenberg矩阵A进行QR分解 mqrdecomp.m

% 程序7.7--mqrdecomp.m
function [Q, R] = mqrdecomp(A)
% 功能:用Givens变换对上Hessenberg矩阵A进行QR分解
% 输入: n阶上Hessenberg矩阵A,其中A(i+1, i) = 0，i > 2.
% 输出:变换后的.上Hessenberg形矩阵A.
[n, n] = size(A); Q = eye(n);
for i = 1 : n-1
    xi = A(i, i); xk = A(i+1, i);
    if xk ~= 0
        d = sqrt(xi ^ 2 + xk ^ 2);
        c = xi / d; s = xk / d;
        J = [c, s; -s, c];
        A(i:i+1, i:n) = J * A(i:i+1, i:n);
        Q(1:n, i:i+1) = Q(1:n, i:i+1) * J';
    end
end
R = A;

% 测试用例
A = [2 3 5 7 8; 4 2 3 5 9; 0 8 3 6 2; 0 0 7 1 3; 0 0 0 6 9];
[Q, R] = mqrdecomp(A)

7.8 QR算法 meigqrdm.m

% 程序7.8--meigqrdm.m
function [Iter, D] = meigqrdm(A, eps)
% 用途:用基本QR算法求实方阵的全部特征值.
% 输入: n阶实方阵A，控制精度eps(默认是1.e-5)
% 输出:迭代次数Iter, A的全部特征值D
% 调用函数: mhessen.m, mqrdecomp.m, eig-仅用于1,2矩阵
if nargin < 2, eps = 1e-5; end
[n, n] = size(A);
D = zeros(n, 1); i = n; Iter = 0; % 初始化
A = mhessen(A); % 化矩阵A为Hessenberg形
% 用基本QR算法进行迭代
while(1)
    if n <= 2
        la = eig(A(1:n, 1:n)); D(1:n) = la'; break;
    end
    Iter = Iter + 1;
    [Q, R] = mqrdecomp(A); % 对上Hessenberg矩阵做QR分解
    A = R * Q; % 做正交相似变换
    % 下面的程序段是判断是否终止
    for k = n-1 : -1 : 1
        if abs(A(k+1, k)) < eps
            if n - k <= 2
                la = eig(A(k+1:n, k+1:n));
                j = i - n + k + 1; D(j:i) = la';
                i = j - 1; n = k; break;
            end
        end
     end
end

% 测试用例
A = [-1 2 3; 2 -3 5; 3 5 -2];
[Iter, D] = meigqrdm(A)

第8章 常微分方程的数值解法

8.1 改进欧拉公式 meuler.m

% 程序8.1--meuler.m
function [x, y] = meuler(df, xspan, y0, h)
% 用途:改进欧拉公式解常微分方程y' = f(x, y), y(x0) = yO
% 格式: [xy] = meuler(df, a, b, yO, h) df为函数f(x, y)，xspan为求解
% 区间[x0, xn]，y0为初值y(x0)，b为步长，[x y]返回节点和数值解矩阵
x = xspan(1) : h : xspan(2); y(1) = y0;
for n = 1 : (length(x)-1)
    k1 = feval(df, x(n), y(n));
    y(n+1) = y(n) + h * k1;
    k2 = feval(df, x(n+1), y(n+1));
    y(n+1) = y(n) + h * (k1+k2) / 2;
end

% 测试用例
df = @(x,y)x + y - 1;
[x, y] = meuler(df, [0,0.5], 1, 0.1)
y1 = exp(x) - x
y - y1

8.2 4阶经典龙格库塔公式 m4rkutta.m

% 程序8.2--m4rkutta.m
function [x, y] = m4rkutta(df, xspan, y0, h)
% 用途:4阶经典龙格库塔公式解常微分方程y' = f(x, y), y(x0) = y0
% 格式: [x, y] = m4rkutta(df, xspan, y0, h)
% df为函数f(x, y)表达式，xspan为求解区间[x0, xn]，
% y0为初值,h为步长，x返回节点，y返回数值解
x = xspan(1) : h : xspan(2); y(1) = y0;
for n = 1 : (length(x)-1)
    k1 = feval(df, x(n), y(n));
    k2 = feval(df, x(n) + h/2, y(n) + h/2 * k1);
    k3 = feval(df, x(n) + h/2, y(n) + h/2 * k2);
    k4 = feval(df, x(n+1), y(n) + h * k3);
    y(n+1) = y(n) + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6;
end

% 测试用例
df = @(x,y)x + y - 1;
[x, y] = m4rkutta(df, [0 0.5], 1, 0.1)
y1 = exp(x) - x
y - y1

8.3 四阶亚当斯“预报-校正”格式 m4adams.m

% 程序8.3--m4adams.m
function [x, y] = m4adams(df, xspan, y0, h)
% 用途:四阶亚当斯“预报-校正”格式解常微分方程y'=f(x,y), y(x0)=yO
% 格式: [x, y] = m4adams(df, xspan, yO, h)
% df为函数f(x,y)表达式，xspan为求解区间[x0, xn],
% y0为初值，h为步长，x返回节点，y返回数值解
x = xspan(1) : h : xspan(2);
[x1, y] = m4rkutta(df, [x(1), x(4)], y0, h);
for n = 4 : (length(x)-1)
    p = y(n) + h/24 * (55 * feval(df,x(n),y(n)) - 59 * feval(df,x(n-1),y(n-1)) ...
        + 37 * feval(df,x(n-2),y(n-2)) - 9 * feval(df,x(n-3),y(n-3)));
    y(n+1) = y(n) + h/24 * (feval(df,x(n-2),y(n-2)) - 5 * feval(df,x(n-1),y(n-1))...
        + 19 * feval(df,x(n),y(n)) + 9 * feval(df,x(n+1),p));
end

% 测试用例
df = @(x,y)x + y - 1;
[x, y] = m4rkutta(df, [0 0.5], 1, 0.1);
y1 = sqrt(1 + 2*x);
[x', y', y1']
y - y1

8.4 四阶经典龙格-库塔公式 m4rkodes.m

% 程序8.4--m4rkodes.m
function [x, y] = m4rkodes(df, xspan, y0, h)
% 用途:四阶经典龙格-库塔公式解常微分方程组y'=f(x,y), y(x0)=yO
% 格式: [x, y] = m4rkodes(df, xspan, y0, h)
% df为向量函数f(x,y)表达式，xspan为求解区间[x0,xn]，
% y0为初值向量，h为步长，x返回节点，y返回数值解向量
x = xspan(1) : h : xspan(2);
y = zeros(length(y0), length(x));
y(:, 1) = y0(:);
for n = 1 : (length(x)-1)
    k1 = feval(df, x(n), y(:,n));
    k2 = feval(df, x(n) + h/2, y(:,n) + h/2 * k1);
    k3 = feval(df, x(n) + h/2, y(:,n) + h/2 * k2);
    k4 = feval(df, x(n+1), y(:,n) + h * k3);
    y(:, n+1) = y(:,n) + h*(k1 + 2*k2 + 2*k3 + k4) / 6;
end

% 测试用例
df = @(x,y)[-0.01 * y(1) - 99.99 * y(2); -100 * y(2)];
[x, y] = m4rkodes(df, [0 500], [2 1], 0.02);
[x, y] = m4rkodes(df, [0 500], [2 1], 0.02);
plot(x, y, 'k');
axis([-50 500 -0.5 2]);
text(120, 0.4, 'y(x)');
text(70, 0.1, 'z(x)');

第9章 蒙特卡洛方法简介

9.1 用蒙特卡洛方法计算pi的近似值 montcpi.m

% 程序9.1--montcpi .m
function [api] = montcpi(n)
% 用途:用蒙特卡洛方法计算pi的近似值
format long;
m = 0; X = 2 * rand(n, 2) - 1;% 产生服从均匀分布的随机数
for(i = 1 : n) % 判定是否落于圆内
    if (X(i,1) ^ 2 + X(i,2) ^ 2 <= 1)
        m = m + 1;
    end
end
api = 4 * m / n;

% 测试用例
api = montcpi(500000)
api = montcpi(500000)
api = montcpi(500000)

9.2 用蒙特卡洛方法求解非线性规划问题 montcnlp.m

% 程序9.2--montcnlp.m
function [sol, zstar] = montcnlp(n)
% 用蒙特卡洛方法求解非线性规划问题
a = 0; b = 10; % 试验点下界和上界
% n=1000000; % 试验点个数
r1 = unifrnd(a, b, n, 1);
r2 = unifrnd(a, b, n, 1);
r3 = unifrnd(a, b, n, 1);
sol = [r1(1), r2(1), r3(1)];
zstar = inf;
for i = 1 : n
    x = [r1(i), r2(i), r3(i)];
    nlpc = nlpconst(x);
    if nlpc == 1
        z = mynlp(x);
        if z <= zstar
        	zstar = z; sol = x;
        end
    end
end
% 目标函数
function z = mynlp(x)
z = 1000 - x(1)^2 - 2*x(2)^2 - x(3)^2 - x(1)*x(2) - x(1)*x(3);
% 約束凾数
function nlpc = nlpconst(x)
t1 = 8*x(1) + 14*x(2) + 7*x(3) - 56;
t2 = x(1)^2 + x(2)^3 + x(3)^3 - 25;
if (abs(t1) <= 0.2 && abs(t2) <= 0.2)
    nlpc = 1;
else
    nlpc = 0;
end

% 测试用例
[sol, zstar] = montcnlp(1000000)
[sol, zstar] = montcnlp(1000000)
[sol, zstar] = montcnlp(1000000)

9.3 用蒙特卡洛方法计算定积分 montcint1.m

% 程序9.3--montcint1.m
function s = montcint1(n)
% 用蒙特卡洛方法计算定积分
% n = 100000; % 试验点个数
format long;
a = 0; b = 1; r = rand(n, 1);
f = @(x)4 ./ (1 + x .^ 2);
x = a + (b - a) * r;
y = f(x);
s = (b - a) * sum(y) / n;

% 测试用例
s = montcint1(100000)
s = montcint1(100000)
s = montcint1(100000)

9.4 测试蒙特卡洛方法计算定积分的收敛速度 montc_rate.m

% 程序9.4--montc_rate.m
function montc_rate()
% 测试蒙特卡洛方法计算定积分的收敛速度
% n = 100000; % 试验点个数
N = [50 100 200 500 1000 2000 5000 10000 20000 50000 100000 200000 500000];
a = 0; b = 1; In = 0;
f = @(x)4 ./ (1 + x .^ 2);
I = quad('4./(1+x.^2)', 0, 1, 1.e-16);
for i = 1:length(N)
    r = rand(N(i), 1);
    x = a + (b-a) * r; z = f(x);
    In = (b-a) * sum(z) / N(i);
    y(i) = abs(In-I) ;
end
N = log(N); y = log(y);
plot(N, y, 'k.-')
xlabel('ln N'); ylabel('ln Err');

% 测试用例
montc_rate()

9.5 用蒙特卡洛方法计算高维积分 montcintn1.m

% 程序9.5--montcintn1.m
function In = montcintn1( )
% 用蒙特卡洛方法计算高维积分
format long;
f = @(x)exp(prod(x, 2)); % 定义被积函数表达式，prod(x, 2)表示x的元素按行连乘积
n = 10000; % 选取随即点的个数
x = rand(n, 4); % 随机生成n个四维单位超立方体内的点
In = sum(f(x)) / n;

% 测试用例
In = montcintn1()
In = montcintn1()
In = montcintn1()

9.6 用蒙特卡洛方法计算不规则区域上的高维积分 montcintn2.m

% 程序9.6--montcintn2.m
function In = montcintn2()
% 用蒙特卡洛方法计算不规则区域上的高维积分
format long;
f = @(x)prod(x); % 定义被积函数表达式，prod(x)表示x的元素按列连乘积
n = 100000;  % 迭取随即点的个数
x1 = unifrnd(1, 2, 1, n); % 随机生成区间[1,2]上n个均匀分布随机数
x2 = unifrnd(1, 4, 1, n); % 随机生成区间[1,4]上n个均匀分布随机数
x3 = unifrnd(1, 16, 1, n); % 随机生成区间[1,16]上n个均匀分布随机数
ind = (x2 >= x1) & (x2 <= 2*x1) & (x3 >= x1.*x2) & (x3 <= 2*x1.*x2);
x = [x1; x2; x3];
In = (4-1) * (16-1) * sum(f(x(:,ind))) / n;

% 测试用例
In = montcintn2()
In = montcintn2()
In = montcintn2()