在进行解方程的时候，如下所示方程
其中，相应的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 牛顿迭代法代码讲解，期末突击干货｜数值计算｜牛顿迭代法｜解题步骤，牛顿求根法
https://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算法的学习
