自学内容网 自学内容网
自学内容网 > 正文

牛顿法、L-M算法

🕗 发布于 2024-10-08 14:21 算法  

在进行解方程的时候，如下所示方程

其中，相应的k11、k12、k21、k22都是已知常量，可以见到其是一个非线性方程。关于非线程方程的求解，我看到网上有两种方法，牛顿法与L-M算法。

1.牛顿法

之前貌似学过，学过也忘记了，没有一点点印象

1.单参数

单参数的牛顿迭代法是非常简单的，c++代码一看就会。

2.双参数牛顿迭代法求解方式

单参数的牛顿法是比较容易求取的，但是双参数的牛顿法在网上很少见到，通过查阅资料可以得到，

单参数：

双参数：

我查了一下资料，上面的双参数的写法（chargtp）是错误的，首先是需要将两个方程变为一个方程。

第4.3讲-牛顿法(上)_哔哩哔哩_bilibili第4.3讲-牛顿法(上), 视频播放量 20003、弹幕量 23、点赞数 265、投硬币枚数 143、收藏人数 360、转发人数 101, 视频作者 MathSu, 作者简介 痴迷于数学研究与教学工作的福哥，带你一起探索高等数学、线性代数及最优化方法等课程中的奥秘，期待每一个聪慧的你加盟福哥学习部落！，相关视频：第4.2讲-最速下降法(上)，第4.3讲-牛顿法(下)，第3.5讲-牛顿法，BFGS牛顿法和拟牛顿法，6.3牛顿迭代法，牛顿法/拟牛顿法/BFGS/DFP算法数学推导(最优化理论)，【数值分析】速成牛顿迭代法｜考试宝典｜一定能听懂！！，9 牛顿迭代法代码讲解，期末突击干货｜数值计算｜牛顿迭代法｜解题步骤，牛顿求根法icon-default.png?t=O83Ahttps://www.bilibili.com/video/BV1JT4y1c7wS/?spm_id_from=333.337.search-card.all.click&vd_source=dd69cecd49e74d6c7b057643f172e2c1

感觉这个视频讲的比较好，因此进行一下记录：

依据上述方法写出自己想要的c++代码

错误版：

首先对于上述公式进行整理，得到：

我的大致的牛顿法函数如下所示：

#include <cmath>
#include <iostream>
#include <Eigen/Dense>
#define M_PI 3.14159265358979323846
#include <iostream>
#include <cmath>
#include <Eigen/Dense>
#include <unsupported/Eigen/NonLinearOptimization>

using namespace Eigen;
using namespace std;


// 定义牛顿法求解方程组的函数,前三个参数分别为输入的参数，第四个为迭代次数，第五个位判断条件
void newton_method(double x0, double x1, double x2, int max_iterations,double epsilon,double x_init,double y_init) {

double x = x_init;
double y = y_init;

for (int i = 0; i < max_iterations; ++i) {
//输入方程

double fxy = sin(x) * x0 - cos(y) * x1 - sin(y) * x2;

//第一次求导：
double dfx_dx1 = cos(x) * x0;//给x求导
double dfy_dy1 = sin(y) * x1 - cos(y) * x2;//给y求导

//将第一次求导变为矩阵
Eigen::Matrix<long double, 2, 1> matrix_dfxy1;
matrix_dfxy1 << dfx_dx1,
dfy_dy1;


//第二次求导：下面方程是一个正定矩阵
double dfx_dx2 = -sin(x) * x0;//给x求导
double dfy_dy2 = cos(y) * x1 + sin(y) * x2;//给y求导

//将第二次求导变为正定矩阵
Eigen::Matrix<long double, 2, 2> matrix_dfxy2;
matrix_dfxy2 << dfx_dx2, 0,
0, dfy_dy2;

//构造牛顿方向
Eigen::Matrix<long double, 2, 1> matrix_xy;
matrix_xy = -matrix_dfxy2.inverse() * matrix_dfxy1;

double  delta_x = matrix_xy(0, 0);
double  delta_y = matrix_xy(0, 1);
;

if (abs(delta_x * delta_x + delta_y * delta_y) < epsilon) {
std::cout << "Converged after " << i + 1 << " iterations." << std::endl;
std::cout << "x = " << delta_x << ", y = " << delta_y << std::endl;
return;
}

x = delta_x;
y = delta_y;
}

std::cout << "Failed to converge after " << max_iterations << " iterations." << std::endl;
}

但是这里要进行初值的赋予，我的想法是针对已经转换好的

进行整理：

先画出b-a的图，这个地方需要用到matlab

我的matlab代码：

然后再进行推导

可得

matlab代码

% 设置 a1 的取值范围
a1 = 0:0.01:2*pi;

% 计算每个 a1 对应的 b1 值
b1 = acos((cos(a1) * 1.99362 - sin(a1) * (-0.00171141)) / 1.99334);

% 绘制函数图形
plot(b1, a1);
xlabel('b1');
ylabel('a1');
title('b1 = arccos((cos(a1) * 1.99362 - sin(a1) * (-0.00171141)) / 1.99262)');
grid on;

然后我们对于相应的代码进行两幅图像整合到一起

相应的matlab代码如下所示：

b = 0:0.001:2*pi;
a = asin((cos(b) * 0.00171141 + sin(b) * 1.99362)/ 1.9939);

plot(b, a);
hold on;

a1 = 0:0.001:2*pi;
b1 = acos((cos(a1) * 1.99362 - sin(a1) * (-0.00171141)) / 1.99334);

plot(b1, a1);

xlabel('b / b1');
ylabel('a / a1');
title('Intersection of asin and arccos Functions');
grid on;

% 找到交点的坐标
intersection_x = [];
intersection_y = [];

for i = 1:numel(b)
    for j = 1:numel(b1)
        if abs(b(i) - b1(j)) < 0.01 && abs(a(i) - a1(j)) < 0.01
            intersection_x = [intersection_x, b(i)];
            intersection_y = [intersection_y, a(i)];
        end
    end
end

% 在交点处标记
plot(intersection_x, intersection_y, 'ro');

% 输出交点坐标
intersection_coordinates = [intersection_x; intersection_y];
disp('Intersection Coordinates:');
disp(intersection_coordinates);

可以得到相应的初值：

注意：这个地方可以看到符合条件的点数是非常非常多的，为什么会出现这个现象，原因是判断条件没有设置好

重新调整阈值
 

b = 0:0.01:2*pi;
a = asin((cos(b) * 0.00171141 + sin(b) * 1.99362)/ 1.9939);

plot(b, a);
hold on;

a1 = 0:0.01:2*pi;
b1 = acos((cos(a1) * 1.99362 - sin(a1) * (-0.00171141)) / 1.99334);

plot(b1, a1);

xlabel('b / b1');
ylabel('a / a1');
title('Intersection of asin and arccos Functions');
grid on;

% 找到交点的坐标
intersection_x = [];
intersection_y = [];

for i = 1:numel(b)
    for j = 1:numel(b1)
        if abs(sqrt((b(i) - b1(j))^2+(a(i) - a1(j))^2 ))< 0.0009
            intersection_x = [intersection_x, b(i)];
            intersection_y = [intersection_y, a(i)];
        end
    end
end

% 在交点处标记
plot(intersection_x, intersection_y, 'ro');

% 输出交点坐标
intersection_coordinates = [intersection_x; intersection_y];
disp('Intersection Coordinates:');
disp(intersection_coordinates);

使用牛顿法求解，得到，明显错误。

正确版：

上面错误的原因是我将两个方程强行变成一个方程进行求解，因此出错。当时想到了这个问题，但是看了看网上很少写到两个方程组的。

当出现多元多次方程的时候，不能够将其变为一个方程进行处理，这种处理方式是错误的。应当将每一个方程看成是一个单独的量进行处理。

可以见到这里是非常的慢的，所以继续我们开始第二个L-M算法的学习

2.L-M算法


原文地址：https://blog.csdn.net/m0_47489229/article/details/135569548

免责声明：本站文章内容转载自网络资源，如本站内容侵犯了原著者的合法权益，可联系本站删除。更多内容请关注自学内容网（zxcms.com）！

相关文章

自学内容网
Copyright © 2015-2024