数值方法：非线性函数求根和求函数极值
数值方法：非线性函数求根和求函数极值
非线性函数求根是数学和工程领域中一个常见的问题。求解非线性方程的根，特别是对于那些没有显式解或者难以通过代数手段求解的复杂方程，通常需要采用数值方法。以下是一些常用的非线性函数求根方法以及它们在工程问题中的应用。
常用方法
- 二分法 ：这是一种简单且直观的方法，适用于在闭区间上连续且在该区间的端点取值异号的函数。二分法的基本思想是将区间不断二分，通过判断函数在区间中点的值来缩小根的所在范围。这种方法简单易行，但收敛速度相对较慢。
- 牛顿迭代法 ：这是一种高效的求根方法，它利用函数的导数和函数值来构造一个迭代公式，通过不断迭代逼近函数的根。牛顿迭代法收敛速度快，但需要函数具有连续的导数，且初始值的选取对收敛性有很大影响。
- 弦截法 ：弦截法是牛顿迭代法的一种改进，它用差商代替导数，从而避免了求导的复杂性。弦截法的收敛速度与牛顿迭代法相近，但计算量稍大。
- 抛物线法 ：抛物线法是一种利用二次插值多项式来逼近函数的方法。它通过构造一个通过三个点的抛物线来逼近函数的曲线，并求解抛物线的根作为函数的近似根。这种方法在某些情况下比二分法和迭代法更高效。
工程应用
非线性函数求根方法在工程问题中有着广泛的应用。例如，在结构力学中，经常需要求解结构的平衡方程，这些方程往往是非线性的。通过求解这些方程的根，可以得到结构的稳定状态或失稳条件。此外，在电路分析、流体力学、热力学等领域中，也经常需要求解非线性方程来得到某些物理量的值。
以电路分析为例，当电路中包含非线性元件（如二极管、晶体管等）时，电路的电压和电流关系将呈现非线性特性。为了分析电路的工作状态，需要求解非线性电路方程。这时，可以利用上述的非线性函数求根方法来求解电路的电压和电流值。
总之，非线性函数求根方法在数学和工程领域中具有广泛的应用价值。掌握这些方法对于解决实际问题具有重要意义。
为了展示非线性方程求根的复杂性和牛顿法等方法的应用，我们可以考虑一个保险公司的利润模型。假设保险公司要确定一个最优的保费(P)，以最大化其长期利润。这个利润不仅依赖于当前的保费和销售量，还依赖于市场份额、客户忠诚度、竞争对手的策略等多个因素。为了简化，我们考虑一个包含指数项和对数项的非线性利润函数。求解的未知数通常是保费水平，即我们希望找到一个保费值，这个值能够使得保险公司在考虑各种风险因素后既能盈利又能吸引客户。
利润函数可以定义为：
%20\times%20P%20-%20\text{Cost}(P)%20-%20\text{Risk}(P))
其中，
- 是保费为 P 时的销售量，可以是一个复杂的非线性函数，如 ，其中 (a, b, c, d) 是参数。
- 是与保费相关的运营成本，也可以是一个非线性函数，如 ，其中 (k, l, m) 是参数。
- 是与保费相关的风险成本，可以假设为 ，其中 是参数。
将上述函数代入利润函数中，我们得到：
)%20\times%20P%20-%20(kP%20+%20lP^2%20+%20mP^3)%20-%20n%20\cdot%20e^{oP})
为了找到最优保费，我们需要对利润函数求导并令其等于0：
这会导致一个非常复杂的非线性方程，其中包含指数、对数和多项式项。这个方程的解不能通过简单的代数方法获得，因此需要使用数值方法，如牛顿法。
牛顿法的基本思想是利用函数的导数和二阶导数（海森矩阵）来逼近函数的零点。对于上述的非线性方程，我们可以定义：
%20=%20\frac{d\text{Profit}}{dP})
然后使用牛顿法的迭代公式：
}{f%27(P_k)})
其中， 是第 k 步的迭代值，) 是 ) 在 处的导数。
通过多次迭代，我们可以找到一个近似解 ，使得 %20\approx%200) ，即利润函数在该点取得极值（可能是最大值或最小值）。为了验证这个点是否是最大值，我们还需要检查利润函数的二阶导数（即海森矩阵）。
请注意，这个模型是高度简化和假设性的，仅用于说明非线性方程求根和优化的复杂性。在实际应用中，模型会更加复杂，并且需要考虑更多的实际因素和数据。此外，使用牛顿法等数值方法时也需要考虑初始值的选择、收敛性和稳定性等问题。
下面提供一组示例参数值，并编写一个简单的北太天元（兼容MATLAB）代码来使用牛顿法求解这个非线性方程。
close all
% 参数定义
a = 1000; % 销售量函数中的参数 a
b = 0.1; % 销售量函数中的参数 b
c = 50; % 销售量函数中的参数 c
d = 1; % 销售量函数中的参数 d
k = 10; % 运营成本函数中的参数 k
l = 2; % 运营成本函数中的参数 l
m = 0.1; % 运营成本函数中的参数 m
n = 50; % 风险成本函数中的参数 n
o = 0.05; % 风险成本函数中的参数 o
figure(1)
P=0:0.1:20;
plot(P, profit(P, a,b,c,d,k,l,m,n,o));
xlabel("保费");
ylabel("收益");
title("北太天元研究保费定价");
figure(2)
%optimal_premium = newton_method_for_premium(a,b,c,d,k,l,m,n,o)
% 调用newton_method_for_premium函数，并获取迭代过程中的节点
[optimal_premium, P_iter, f_prime_iter] = newton_method_for_premium_plot(a,b,c,d,k,l,m,n,o);
% 绘制利润曲线
P=0:0.1:20;
plot(P, profit(P, a,b,c,d,k,l,m,n,o));
hold on; % 保持图像，以便在同一图上绘制迭代点
xlabel("保费");
ylabel("收益");
title("北太天元研究保费定价与牛顿法迭代过程");
% 绘制牛顿法迭代过程
plot(P_iter, profit(P_iter, a,b,c,d,k,l,m,n,o), 'ro-'); % 迭代点的收益
plot(P_iter, zeros(size(P_iter)), 'bx-'); % 迭代点在x轴上的投影
legend('利润曲线', '迭代点收益', '迭代过程');
figure(3)
% 假设您已经有了迭代点的坐标 P_iter (x坐标) 和 f_iter (y坐标)
x_iter = P_iter; % 迭代点的x坐标（保费值）
y_iter = profit(P_iter, a,b,c,d,k,l,m,n,o); % 假设 f_profit 是利润函数，计算每次迭代的利润
% 绘制利润曲线（示例，您应该使用您自己的利润函数）
x = linspace(0, max(P_iter), 100);
y = profit(x, a,b,c,d,k,l,m,n,o); % 利润函数对应的y值
plot(x, y, 'b-', 'LineWidth', 2); % 绘制利润曲线
hold on;
% 绘制迭代点并标注序号
for i = 1:length(P_iter)
plot(x_iter(i), y_iter(i), 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r'); % 绘制点
% 在点附近添加序号文本
text(x_iter(i), y_iter(i), sprintf('{%d}', i), ...
'VerticalAlignment', 'bottom', ...
'HorizontalAlignment', 'right', ...
'FontSize', 8);
pause(1)
end
% 设置图表标题和坐标轴标签
title('Newton Method Iteration with Iteration Numbers');
xlabel('Premium (P)');
ylabel('Profit');
grid on;
hold off;
function [optimal_premium, P_iter, f_prime_iter] = newton_method_for_premium_plot(a,b,c,d,k,l,m,n,o)
% 初始保费猜测值
P = 10;
% 容忍误差和最大迭代次数
tol = 1e-6;
max_iter = 100;
% 初始化迭代过程记录数组
P_iter = zeros(1, max_iter);
f_prime_iter = zeros(1, max_iter);
% 牛顿法迭代
for iter = 1:max_iter
% 计算一阶导数 f'(P)
f_prime = derivative_profit(P, a, b, c, d, k, l, m, n, o);
% 记录迭代节点
P_iter(iter) = P;
f_prime_iter(iter) = f_prime;
% 检查是否达到收敛条件
if abs(f_prime) < tol
break;
end
% 计算二阶导数 f''(P)
f_double_prime = second_derivative_profit(P, a, b, c, d, k, l, m, n, o);
% 检查二阶导数是否为零，以避免除以零的错误
if abs(f_double_prime) < tol
disp('Second derivative is too small, stopping iteration.');
break;
end
% 更新保费值
P = P - f_prime / f_double_prime;
end
% 截断未使用的迭代记录
P_iter = P_iter(1:iter);
f_prime_iter = f_prime_iter(1:iter);
% 输出最优保费
optimal_premium = P;
end
function optimal_premium = newton_method_for_premium(a,b,c,d,k,l,m,n,o)
% 初始参数设置
%a = 1000; b = 0.1; c = 50; d = 1; k = 10; l = 2; m = 0.1; n = 50; o = 0.05;
% 初始保费猜测值
P = 100;
% 容忍误差和最大迭代次数
tol = 1e-6;
max_iter = 100;
% 牛顿法迭代
for iter = 1:max_iter
% 计算一阶导数 f'(P)
f_prime = derivative_profit(P, a, b, c, d, k, l, m, n, o);
% 计算二阶导数 f''(P)
f_double_prime = second_derivative_profit(P, a, b, c, d, k, l, m, n, o);
% 检查二阶导数是否为零，以避免除以零的错误
if abs(f_double_prime) < tol
disp('Second derivative is too small, stopping iteration.');
break;
end
% 更新保费值
P = P - f_prime / f_double_prime;
% 检查收敛性：当一阶导数足够接近0时，我们认为找到了极值点
if abs(f_prime) < tol
break;
end
end
% 输出最优保费
optimal_premium = P;
end
function f = profit(P, a, b, c, d, k, l, m, n, o)
% 销售量函数
N = a * exp(-b * P) + c * log(P + d);
% 运营成本函数
Cost = k * P + l * P.^2 + m * P.^3;
% 风险成本函数
Risk = n * exp(o * P);
% 利润函数
f = N .* P - Cost - Risk;
end
function df = derivative_profit(P, a, b, c, d, k, l, m, n, o)
% 利润函数的一阶导数
dN_dP = -a * b * exp(-b * P) + c / (P + d);
df = (a * exp(-b * P) + c * log(P + d)) + P * dN_dP - (k + 2 * l * P + 3 * m * P^2) - n * o * exp(o * P);
end
function ddf = second_derivative_profit(P, a, b, c, d, k, l, m, n, o)
% 利润函数的二阶导数
ddf = -2*a*b*exp(-b*P) + c/(P + d) + a*b^2*P*exp(-b*P) + c*d/((P + d)^2) - 2*l - 6*m*P - n*o^2*exp(o*P);
end
收益作为保费的函数图像,从中可以看出最大收益的保费在10到14之间.
newton法迭代,得到收益最大的保费是12.2731, 迭代中的点也画在图中
这是动画的最后一帧,每隔一秒画一个迭代点