解非线性方程组
实验类型：●验证性实验 ○综合性实验 ○设计性实验
实验目的：进一步熟练掌握解非线性方程组牛顿迭代算法，提高编程能力和解算非线性方程组问题的实践技能。
实验内容：
设有非线性方程组(此方程组是非标准型)
实验说明：要求先严格按照牛顿迭代法解非线性方程组的算法步骤（参照课堂课件PPT）手工计算非线性方程组的解—手工迭代计算前三步，把计算及循环结束判定等详细过程写在A4规格的打印纸上；画出牛顿迭代法解非线性方程组的算法程序设计流程图；编写程序并用计算机计算、输出中间结果和最终结果，列表比较计算机计算中间结果（前三个近似解向量）与手工计算（前三个近似解向量）是否一致。列表对比后，给出比较后的验证结论。
实验步骤
- 要求上机实验前严格按照牛顿迭代法解非线性方程组的算法步骤（参照课堂PPT）对于非线性方程组标准化成F(P)=(f1(P), f2(P))T=(f1(x1, x2), f2(x1, x2))T=(x1,x2) T =0手工迭代计算三步（可使用计算器），并把计算过程及其循环结束判定过程步骤详细写在A4纸上。
- ) T =0手工迭代计算三步（可使用计算器），并把计算过程及其循环结束判定过程步骤详细写在A4纸上。
- 对手工计算做进一步的总结梳理，然后画出牛顿迭代法解非线性方程组的算法程序设计流程图；流程图必须符合流程图的图符规范。应用Matlab语言在纸上编写出牛顿迭代法解非线性方程组算法的函数程序代码。
- 编辑录入函数程序。
- 输入数据（截图），调试程序并记录调试过程中出现的问题（截图）及修改程序的过程（截图），把你的截图及时粘贴在Word文档中，并对截图配上必要的文字说明，以便证明和描述你的实验过程真实可靠。
- 经反复调试后，运行程序并验证程序运行的前三步中间结果和手工计算的前三步中间结果是否一致，要求通过列表对比、分析和论证程序运行的前三步中间结果和手工计算的前三步中间结果是否一致，若不一致，必须找到不一致的原因，并预以纠正，最后明确给出对比分析结论 。
- 将手工计算过程及其验证结论拍照成JPG图片插入在实验报告的word 电子版中。
数学建模阅读与思考：利用利用校园网-陕西科技大学-图书馆-国内数据库（http://library.sust.edu.cn/gnskj/index.jhtml），自己在“cnki知识网络服务平台”或者“万方数据库”查阅附件给定的数学建模文献资料、阅读文献资料并写出分析建模的详细过程、求解计算过程、实际问题解答或阐释的内容。了解牛顿迭代法解非线性方程组的应用。
实验报告书写要求：根据实验情况和结果撰写并递交实验报告。实验报告应当有算法原理简介，算法流程图，程序代码，运行调试记录，计算结果阐释；对于数学建模的应用问题要求利用校园网-陕西科技大学-图书馆-国内数据库（http://library.sust.edu.cn/gnskj/index.jhtml），自己在“cnki知识网络服务平台”或者“万方数据库”查阅附件给定的文献资料、阅读文献资料并写出分析建模的详细过程、求解计算过程、实际问题解答或阐释的内容。在A4纸上对给定的非线性方程组和文献上的非线性方程组问题用牛顿迭代法手工迭代计算前三步，手工计算要求所有数据保留6位小数，然后将手工迭代计算前三步的计算结果与计算机程序计算的前三步结果进行比较，若数据整数及小数点后4位小数都相同则认为是相同的，若不同，要分析不同的原因并予以校正，经分析、校正后给出比较的结论。 将手工计算过程及其验证结论拍照成JPG图片插入在实验报告的word 电子版中。
实验报告打印和装订顺序要求
实验报告打印要求：在A4纸上将实验报告的Word文档双面打印；打印的字迹要清晰，报告内容文字或截图上的文字打印后，要确保字迹清晰可见，若字迹不清影响评阅，报告得分将酌情扣减。
实验报告装订顺序
- 实验任务（本文档版式不得做任何改动）；
- 实验报告正文（Word文档）；
- 可以将A4纸上详细手工计算过程拍照成JPG图片插入在实验报告的word 电子版中随Word文档双面打印。
- 左边沿装订
实验总结(学会了......; 掌握了......; 训练了......; 发现了......; 今后学习中......有待提高。）
程序代码
电子报告word文件命名规则：专业班级-完整学号-实验X-姓名.doc, 如信息123班学号为201212030315的郭海涛同学实验2报告word文件命名则应是：信息123-201212030315-实验2-郭海涛.docx, 其中 .docx是Word文件扩展名。
特别提醒：学委负责检查电子文件命名规范，命名不规的将不接收。每人必需提交电子报告和纸质报告各一份(两者内容一致)。本实验任务书版式不得做任何改动。电子报告以班为单位提交压缩包.
附件：数学建模参考文献材料：
朱喜明，王存良*
（中国电子科技集团第27研究所，郑州 450005）
摘 要：给出了用户接收机只收到三颗导航卫星信号时的二维定位解算模型及 其解算方法，重点介绍了解算方法中接收机概略坐标的求法及利用得到的概略坐标 求解接收机位置和钟差，并做了仿真解算。从卫星和接收机不同位置的很多次仿真 解算结果表明，接收机位置和钟差能够极收收敛，通常迭代两次就能得到非常好的结 果，完全适合工程应用。
关键词：二维定位；概略坐标；迭代求解
| 此页及其前的内容为实验任务内容，实验任务版式不得做任何更改；实验报告从下页开始编写排版。 |
算法原理：
手工迭代计算：
（1）x0 = [2.00; 0.25];
（2）x0 = [-0.20; 1.00];
（3）x0 = [-1.00; -1.00];
程序代码：
function main()
% 初始值
P0 = [2.00; 0.25];
fprintf('初始迭代值(2.00, 0.25):\n');
newton_method(P0);
P0 = [-0.20; 1.00];
fprintf('初始迭代值(-0.20, 1.00):\n');
newton_method(P0);
P0 = [-1.00; -1.00];
fprintf('初始迭代值 (-1.00, -1.00):\n');
newton_method(P0);
end
function newton_method(P0)
F = @(P) [P(1)^2 - 2*P(1) - P(2) + 0.5; P(1)^2 + 4*P(2)^2 - 4];
J = @(P) [2*P(1)-2, -1; 2*P(1), 8*P(2)];
tol = 1e-6;
max_iter = 100;
P = P0;
for i = 1:max_iter
% 计算函数值和雅可比矩阵
f = F(P);
Jf = J(P);
% 计算更新量
delta_P = -inv(Jf) * f;
% 更新变量
P = P + delta_P;
% 判断收敛
if norm(delta_P) < tol
fprintf('经过 %d 次迭代\n', i);
fprintf('最后: x1 = %.6f, x2 = %.6f\n', P(1), P(2));
break;
end
fprintf('第%d次迭代近似值为: x1 = %.6f, x2 = %.6f\n',i, P(1), P(2));
end
% 判断发散
if i == max_iter
fprintf("牛顿迭代法发散\n");
end
end
手工计算结果与程序运行结果前三次比对分析：
| （1）x0=[2.00;0.25]; | ||||
| 迭代次数 | 手工x1 | Matlab的x1 | 手工x2 | Matlab的x2 |
| 第一次 | 1.906250 | 1.906250 | 0.312500 | 0.312500 |
| 第二次 | 1.900691 | 1.900691 | 0.311213 | 0.311213 |
| 第三次 | 1.900677 | 1.900677 | 0.311219 | 0.311219 |
| （2）x0=[-0.20;1.00]; | ||||
| 迭代次数 | 手工x1 | Matlab的x1 | 手工x2 | Matlab的x2 |
| 第一次 | -0.222449 | -0.222449 | 0.993878 | 0.993878 |
| 第二次 | -0.222215 | -0.222215 | 0.993808 | 0.993808 |
| 第三次 | -0.222215 | -0.222215 | 0.993808 | 0.993808 |
| （3）x0=[-1.00;-1.00]; | ||||
| 迭代次数 | 手工x1 | Matlab的x1 | 手工x2 | Matlab的x2 |
| 第一次 | 0.166667 | 0.166667 | -1.166667 | -1.166667 |
| 第二次 | 0.873543 | 0.873543 | -0.983683 | -0.983683 |
| 第三次 | 1.756130 | 1.756130 | -0.707227 | -0.707227 |
我将手工计算和程序计算的前三步结果进行了比较，发现整数及小数点后4位小数都相同，说明程序计算的结果与手工计算结果一致。
数学建模：
1.分析建模
给出了用户接收机只收到三颗导航卫星信号时的二维定位解算模型及其解算方法，重点介绍了解算方法中接收机概略坐标的求法及利用得到的概略坐标求解接收机位置和钟差，并做了仿真解算。从卫星和接收机不同位置的很多次仿真
解算结果表明，接收机位置和钟差能够极收收敛，通常迭代两次就能得到非常好的结果，完全适合工程应用。
为了简化计算，假定已利用了导航电文中的卫星星历表预报参数，算出了三颗卫星在地心空间直角坐标系中的坐标分别为（x1，y1，z1）、（x2，y2，z2）、 （x3，y3，z3）海平面上某一位置在地心空间直角坐标系中的坐标为（x，y，z），设在该位置的用户接收机观测三颗星得到伪距观测值，且已进行了各项误差修正（包括卫星钟钟差改正），接收机时钟相对该导航系统时间基准的钟差为ΔtR；则可得观测方程
ρj＝Rj＋cΔtR（j＝1,2,3）
可写为
ρj＝［（x－xj）2＋（y－yj）2 +（z－zj）2］1／2＋d…(1)
式中d＝cΔtR,为接收机钟差的等效距离偏差。
由于接收机在海面上,所以其坐标满足方程
式中a＝6378140m,b＝6356755．28816m
解（1）、（2）组成的方程组,可得x,y,z,d,进而算出接收机所在位置的经度L、纬度B和钟差ΔtR。经仿真解算证明该方程组正确可行,满足实用需求。（如果x,y,z用经纬度表示出来,方程组（1）就变成了关于L、B、d的超越方程组,非常难解）
2. 求解计算过程
上述定位方程组的解法,采用迭代求解。迭代求解首先要解决迭代初始解问题,这里由于d的符号不能事先确定,所以采用先求出接收机的概略坐标(接收机概略坐标指的是接收机在地心空间直角坐标系中的大概位置坐标),并作为初值求出概略坐标的修正量和d;将概略坐标加上其修正量作为第二次的迭代值再重新解算,直到概略坐标的修正量满足要求为止。
2)计算
3）计算
A＝m22＋n22＋1，B＝2（m2x1＋n2y1-z1－m1m2－n1n2）
C＝m12＋n12－2m1x1－2n1y1-ρ12＋x12＋y12+z12
若B2-4AC<0，则方程组（1）无实数解，终止求解过程。此种情况无法定位，该导航系统通常会通过其完好性监测告知用户。
地球表面上二维定位仿真计算。由题易知三颗卫星的高度皆为20200km，某一时刻其经纬度分别为40°、112°，35°、114°，32°、117°;用户接收机在海平面上某一位置的经纬度假定为35°、115°。之所以先假定出接收机的位置 ，是为了与仿真计算得到的接收机位置作比较，观察算法误差。
通过MATLAB牛顿迭代法计算，三颗卫星的地心大地坐标转换为地心空间直角坐标分别为：
x1=-7629540.435m,y1=18883775.23m,z1=17062297.19m;
x2=-8857626.141m,y2=19894554.04m,z2=15224112.62m;
x3=-10235048.89m,y3=20087414.48m,z3=14064802.14m。
将x1、y1、z1的值代入方程所得的值非常接近1;将x2、y2、z2的值代入方程的左端所得的值远大于1;所以x1、y1、z1的的值是接收机的概略坐标。
3．实际问题解答或阐释的内容
经过大量仿真表明：对于大地高在45m以下的地区，可以按照海平面上的二维定位进行计算，其定位误差随着大地高的减小而变小；对于大地高在45m以上的地区，按照海平面上的二维定位进行计算，通常定位误差比概略位置误差还大，所以不再应用地球椭球方程，应采用方程组进行定位，将得到的两组解分别转换为地心大地坐标系中的坐标，其大地高的绝对值较小的那组解即为接收机所在的位置，由于忽略了接收机的钟差，这种情况定位误差较大。
实验总结：
学会了使用牛顿迭代法求解非线性方程组的基本原理和方法。通过实验，我了解到牛顿迭代法是一种有效的数值计算方法，可以用于求解复杂的非线性方程组。
掌握了编写 MATLAB 程序实现牛顿迭代法的能力。在实验中，我学会了如何利用 MATLAB 编写程序来求解非线性方程组，提高了我的编程能力。
训练了分析和解决实际问题的能力。实验要求结合校园网数据库中的文献资料进行数学建模分析，这锻炼了我的问题分析和求解能力，提高了我的实际应用能力。
发现了牛顿迭代法的收敛性和计算精度的影响因素。在实验中，我发现初始值的选择对于牛顿迭代法的收敛性和计算精度有很大影响，这启示我在实际应用中要注意初始值的选择。
今后学习中需要进一步提高数学建模和算法设计的能力。通过本次实验，我意识到数学建模和算法设计是非常重要的能力，今后我将继续学习和提高这方面的能力，以更好地应用于实际问题的求解中。