matlab：涉及复杂函数图像的交点求解

news2024/10/6 18:35:01

matlab：涉及复杂函数图像的交点求解

在MATLAB中求解两个图像的交点是一个常见的需求。本文将通过一个示例，展示如何求解两个图像的交点，并提供相应的MATLAB代码。

画出图像

首先，我们需要绘制两个图像，以便直观地看到它们的交点。以下是绘制图像的MATLAB代码：

% 定义符号变量
syms x1 x2;

% 上边界方程
eq1 = 10 + 110 * (0.8 + 0.05 + 0.4 * sin(4 * atan2(x2, x1))^16)^2 - (x1 + x2) == 0;

% 绘制图形
figure;

% 使用 fimplicit 绘制上边界
fimplicit(@(x1, x2) 10 + 110 * (b + 0.05 + 0.4 * sin(4 * atan2(x2, x1))^16)^2 - (x1 + x2), [0, 100, 0, 100]);
hold on;

% 使用 fimplicit 绘制 y = 100 - x1
fimplicit(@(x1, x2) x1 + x2 - 100, [0, 100, 0, 100]);

% 设置图例和标题
legend('Upper Boundary', 'y = 100 - x1');
title('Plot of Equations');
xlabel('x1');
ylabel('x2');
hold off;

绘制的图像如下所示：
在这里插入图片描述
从图中可以看出，存在四个实数域交点。

求解交点方法1：solve不加求解范围（失败）

首先尝试使用solve函数来求解交点：

% 定义符号变量
syms x1 x2;

% 定义参数
b = 0.8;
l = atan2(x2, x1);

% 定义方程组
eq1 = 10 + 110 * (b + 0.05 + 0.4 * sin(4 * l)^16)^2 - (x1 + x2);
eq2 = x1 + x2 - 100;

% 解方程组
[sol_x1, sol_x2] = solve([eq1 == 0, eq2 == 0], [x1, x2]);

% 转换为数值解
sol_x1 = double(sol_x1);
sol_x2 = double(sol_x2);

% 筛选实数解
real_solutions = [sol_x1, sol_x2];
real_solutions = real_solutions(imag(real_solutions(:, 1)) == 0 & imag(real_solutions(:, 2)) == 0, :);

% 输出实数解
disp('Real solutions (x1, x2):');
disp(real_solutions);

输出为空值，这表明在解的过程中遇到了问题。

求解交点方法2：solve加求解范围（成功）

% 定义符号变量
syms x1 x2;

% 求解范围
assume(x1>=0&x1<=100)
assume(x2>=0&x2<=100)

% 定义参数
b = 0.8;
l = atan2(x2, x1);

% 定义方程组
eq1 = 10 + 110 * (b + 0.05 + 0.4 * sin(4 * l)^16)^2 - (x1 + x2);
eq2 = x1 + x2 - 100;

% 解方程组
[sol_x1, sol_x2] = solve([eq1 == 0, eq2 == 0], [x1, x2]);

% 转换为数值解
sol_x1 = double(sol_x1);
sol_x2 = double(sol_x2);

% 筛选实数解
real_solutions = [sol_x1, sol_x2];
real_solutions = real_solutions(imag(real_solutions(:, 1)) == 0 & imag(real_solutions(:, 2)) == 0, :);

% 输出实数解
disp('Real solutions (x1, x2):');
disp(real_solutions);

输出：
在这里插入图片描述

求解交点方法3：fsolve（成功）

下面通过数值方法而不是符号方法来找到解，使用 fsolve（数值求解函数）：

% 定义匿名函数
func = @(x) [10 + 110 * (b + 0.05 + 0.4 * sin(4 * atan2(x(2), x(1)))^16)^2 - (x(1) + x(2)), x(1) + x(2) - 100];

% 设置选项以使用较大的初始搜索范围
options = optimoptions('fsolve', 'Display', 'off', 'MaxFunctionEvaluations', 6000, 'MaxIterations', 4000);

% 存储解
solutions = [];

% 尝试多个随机初始猜测
for i = 1:100
    initial_guess = rand(1, 2) * 100; % 生成0到100之间的随机初始猜测
    [sol, fval, exitflag, output] = fsolve(func, initial_guess, options);
    % 只有当fsolve成功收敛时才记录解
    if exitflag > 0 && all(abs(fval) < 1e-6)
        solutions = [solutions; sol];
    end
end

% 去除重复的解，考虑数值误差
solutions = round(solutions, 3); % 四舍五入到三位小数
solutions = unique(solutions, 'rows', 'stable');

% 过滤掉不在感兴趣区域的解
solutions = solutions(all(solutions >= 0 & solutions <= 100, 2), :);

% 输出数值解
disp('Numerical solutions (x1, x2):');
disp(solutions);