当前位置：首页 > 运维知识 > windows >内容正文

windows

线性系统实验：化学方程式配平与天体轨道参数估计

发布时间：2023/12/20 windows 52 豆豆

生活随笔收集整理的这篇文章主要介绍了线性系统实验：化学方程式配平与天体轨道参数估计小编觉得挺不错的,现在分享给大家,帮大家做个参考.

线性系统

消元法
- 高斯消元法
- 高斯-约旦消元法
- - 工程应用：化学方程式配平
逆矩阵求解
- 矩阵的逆
解的结构
- 超定方程组近似估值
- - 工程应用：天体轨道参数估计

消元法

线性系统，在高中称为 多元一次方程组，因为在线性代数里，我们把矩阵看成系统，而这些方程组的未知数都只有一次，所以就了线性系统。

${1x+2y=53x+4y=6\left\{\begin{matrix} 1x + 2y = 5 & \\ 3x + 4y = 6 & \end{matrix}\right.$

我们把真实世界的问题，转换为线性方程（线性系统），研究线性系统，也就是解这些方程，解出来问题就解决了。

高斯消元法

初中的时候，我们学习了高斯消元法。

高斯消元法主要分为两步，

消元：要减少某些方程中元的数量；
回代：是把已知的解代入到方程式中，求出其他未知的解。

如果方程和元的数量很小，那么高斯消元法并不难理解。

可是如果方程和元的数量很多，整个过程就变得比较繁琐了。

而后，自然而然人们想到了对计算进行批处理，可以把高斯消元法转为矩阵的操作。

比如：

${1x+2y+4z=73x+7y+2z=−112x+3y+3z=1\left\{\begin{matrix} 1x + 2y + 4z = 7 ~~~~\\ ~3x + 7y + 2z = -11\\ 2x + 3y + 3z = 1~~~~\\ \end{matrix}\right.$

用矩阵表示(矩阵乘法)：

$[124372233]∗[xyz]=[7−111]\left[ \begin{matrix} 1 & 2 & 4 \\ 3 & 7 & 2 \\ 2 & 3 & 3 \end{matrix} \right]*\left[ \begin{matrix} x\\ y\\ z \end{matrix} \right]=\left[ \begin{matrix} 7\\ -11\\ 1 \end{matrix} \right]$

我们再把系数矩阵和等式右边的矩阵放在一起：

${1247372−112331}\left\{ \begin{array}{ccc|c} 1 & 2 & 4 &7 \\ 3 & 7 & 2 & -11 \\ 2 & 3 & 3 & 1 \end{array} \right\}$

这个矩阵叫【增广矩阵】，就是在系数矩阵右边增加了一列等号右边的解。

通过矩阵解线性方程，本质也是消元法，要做的就是把消元的过程封装到矩阵里。

用消元法解方程，因为操作的是矩阵，所以消元法的书写方式有点变化：

原【一个方程的左右俩边同时乘以一个常数】变为【矩阵的某一行乘以一个常数】
原【一个方程加（或减）另一个方程】变为【矩阵的一行加（或减）另一行】
原【交换俩个方程的位置】变为【交换矩阵的俩行】

也就是说，计算方式是完全一致的，只不过操作对象从方程变成了矩阵，因为矩阵是可以批处理运算的工具。

我快速手算一下：

从上到下，以矩阵的第一行为基准，消去第二、第三行的第一个元素（系数化为 $0$ ）

第一步 $(2) - (1) * 3$ ，把第二行的第一个元素化成 $0$ ；

第二步 $(3) - (1) * 2$ ，把第三行的第一个元素化成 $0$ ；

第三步 $(3) * - 1 - (1)$ ，把第三行的第二个元素化成 $0$ ；

第三步后矩阵的第三行是 $[00−1545]\left[\begin{array}{ccc|c}0 & 0 & -15 & 45\\ \end{array} \right]$ ，而后给第三行左右俩边同时除以 $- 15$ ，最后的结果就是上图最后一个矩阵。

现在我们知道，最后一个未知数 $z = 3$ ，这是消元法的消元部分，接着是回代部分，由下往上回代解出其余未知数：）。

总结一下，高斯消元法。

消元部分：

从上到下，以矩阵的第一行为基准，消去第二、第三行的第一个元素（系数化为 $0$ ），这就代表消去了一个元；

因为他们都是基于矩阵第一行第一个位置消元的，所以这个位置也被称为【主元】。

不仅第一行第一个位置是主元，第二行第二个、第三行第三个、第 $n$ 行第 $n$ 个都是主元。

P.S. 只要在这个位置，任何不为 0，就可以当主元。如果第一行第一个元素为 $0$ ，就需要【交换矩阵的俩行】。若第一行主元不为 $1$ 的元素，先化为 $1$ ，再运算。

首先就是把这些主元位置的系数化为 $1$ ，之后其他行就可以使用【矩阵的某一行乘以一个常数】、【矩阵的一行加（或减）另一行】把每行主元左边的元素化为 $0$ 。

最后，矩阵会变成一个阶梯型的矩阵，比如这样。

回代部分：

由下往上回代解出其余未知数，把最后一个未知数代入倒数第二行，得到倒数第二个未知数的解；再把倒数第一、倒数第二的未知数代入到第一个里面，得到 $3$ 个未知数的解。

高斯消元法实现：

// 运行：命令行输入 gcc/g++ 当前源文件.c/cpp #include<stdio.h> #include<stdlib.h> #include<math.h>double **A = NULL, *b = NULL, *x = NULL; unsigned int RANK = 4;unsigned int makematrix(){unsigned int r, c;printf("请输入矩阵行列数，用空格隔开：");scanf("%d %d", &r, &c);A = (double**)malloc(sizeof(double*)*r);// 创建一个指针数组，把指针数组的地址赋值给a ,*r是乘以r的意思for (int i = 0; i < r; i++)A[i] = (double*)malloc(sizeof(double)*c);// 给第二维分配空间for (int i = 0; i < r; i++) {for (int j = 0; j < c; j++)A[i][j] = 0.0; }b = (double*)malloc(sizeof(double)*r);for (int i = 0; i < r; i++){b[i] = 0.0;}x = (double*)malloc(sizeof(double)*c);for (int i = 0; i < c; i++){x[i] = 0.0;}return r;// 一般都是输入方阵，返回行数也阔以 }void getmatrix(void) { // 输入矩阵并呈现printf("\n\n请按行从左到右依次输入系数矩阵A，不同元素用空格隔开：\n");for (int i = 0; i < RANK; i++) {for(int j = 0;j<RANK;j++) {scanf("%lf", &A[i][j]);}}printf("\n\n系数矩阵如下：\n");for (int i = 0; i < RANK; i++) {for (int j = 0; j<RANK; j++) {printf("%g\t",A[i][j]);}putchar('\n');}printf("\n\n请按从上到下依次输入常数列b，不同元素用空格隔开：\n");for (int i= 0; i<RANK; i++) {scanf("%lf", &b[i]);}printf("常数列如下\n");for (int i = 0; i<RANK; i++) {printf("%g\t", b[i]);}putchar('\n'); }void Gauss_calculation(void) { // Gauss消去法解线性方程组double get_A = 0.0;printf("\n\n利用以上A与b组成的增广阵进行高斯消去法计算方程组\n");for (int i = 1; i < RANK; i++) {for (int j = i; j<RANK; j++) {get_A = A[j][i - 1] / A[i - 1][i - 1];b[j] = b[j] - get_A * b[i - 1];for (int k = i-1; k < RANK; k++) {A[j][k] = A[j][k] - get_A * A[i-1][k];}}}printf("\n\n顺序消元后的上三角系数增广矩阵如下：\n");for (int i = 0; i < RANK; i++) {for (int j = 0; j<RANK; j++) {printf("%g\t", A[i][j]);}printf(" %g", b[i]);putchar('\n');}printf("\n\n利用回代法求解上三角方程组，解得：\n\n");for (int i = 0; i < RANK; i++) {double get_x = 0.0;for (int j = 0; j < RANK; j++) {get_x = get_x + A[RANK-1-i][j]*x[j];// 把左边全部加起来了，下面需要多减去一次Xn*Ann}x[RANK - 1 - i] = (b[RANK - 1 - i] - get_x + A[RANK - 1 - i][RANK - 1 - i] * x[RANK - 1 - i]) / A[RANK - 1 - i][RANK - 1 - i];}for (int i = 0; i < RANK; i++) {printf("x%d = %5g\n", i + 1, x[i]);}printf("\n\n计算完成，按回车退出程序或按1重新输入矩阵\n"); }int main() {RANK = makematrix();getmatrix();Gauss_calculation();// P.S. double **A, *b, *x 没释放return 0; }

高斯消元法演示：

高斯-约旦消元法

高斯消元法的问题，只能通过矩阵得到最后一个未知数的解，之后还得一直回代依次得到所有结果，那有木有什么方法一次性把结果捣鼓出来呢？

有呀，高斯-约旦消元法。

回到高斯消元法的消元后，其实我们也可以不用回代，只有把 $- 10$ 化为 $0$ 即可。

我们倒着进行消元即可，初始的消元是从第一行开始，从上往下以第一行的主元为基准，把主元下面的位置都化为 $0$ ，现在反过来从最后一行开始，从下往上以最后一行的主元为基准，把主元上面的位置都化为 $0$ 。

高斯-约旦消元法：

高斯消元：前向过程，从上到下
约旦消元：后向过程，从下到上

高斯-约旦消元法实现：

// 运行：在命令行输入 g++ -std=c++11 当前源文件.cpp #include<cstdio> #include<cmath> #include<algorithm> #define eps 1e-8using namespace std;double a[55][55], ans[55]; // a为增广矩阵 int d; int gauss_jordan(int n) { int r, w = 0;for (int i = 0; i < n && w < n; w++, i++) { // 进行到第i列，第w行 int r = w;for (int j = w + 1; j < n; j++)if (fabs(a[j][i]) > fabs(a[r][i]))r = j; // 找到当前列绝对值最大的行 if (fabs(a[r][i]) < eps) {w--;continue;} // 当前列值都为0,跨过当前步 if (r != w)for (int j = 0; j <= n; j++)swap(a[r][j], a[w][j]); // 交换当前列绝对值最大的行和没计算过的第一行，使用最大值运算，可减少误差for (int k = 0; k < n; k++) { // 消去当前列（除本行外）if (k != w)for (int j = n; j >= w; j--)a[k][j] -= a[k][i] / a[w][i] * a[w][j];}}return w; }int main() {int n;scanf("%d", &n); // 默认输入的矩阵是方阵，行数等于列数，但这个限制也不是必须的，可以打破for (int i = 0; i < n; i++) {for (int j = 0; j <= n; j++) { // 因为解的那一列也一起输入，所以加了一行scanf("%lf", a[i] + j); // C/C++ 里其实只有一维数组，a[i] + j = a[i][j]}}d = gauss_jordan(n);d--;for (int i = 0; i < n; i++) { // 有一个方程 = 右边不为0 = 左边为0，则无解 bool d = 1;for (int j = i; j < n; j++)d &= (fabs(a[i][j]) < eps);if (d && fabs(a[i][n]) > eps) {putchar('-1');return 0;}}for (int i = 0; i < n; i++) { // 消元后有变量在多个方程中出现，则有多个解 int max1 = 0;for (int j = i; j < n; j++)if (fabs(a[i][j]) > eps)max1++;if (max1 > 1) {putchar('0');return 0;}}for (int i = 0; i < n; i++)ans[i] = a[i][n] / a[i][i];putchar('\n');for (int i = 0; i < n; i++) {if (fabs(ans[i]) < eps)printf("x%d = 0\n", i + 1);elseprintf("x%d = %5.2lf\n", i + 1, ans[i]);}return 0; }

高斯-约旦消元法演示：

工程应用：化学方程式配平

问题描述：给出一个未配平的化学方程式，根据质量守恒定律对其分配，不考虑化合价问题

示例：）

输入：Cu+HNO3=Cu(NO3)2+NO+H2O（中间不要加入空格）
输出：3Cu+8HNO3=3Cu(NO3)2+2NO+4H2O

分析：

化学反应遵循质量守恒定律，反应前后的原子种类和数量保持不变，根据这，采用【待定系数法】来配平化学方程式。

具体的操作，给方程式中的每一项设一个待定系数，列出方程组。

比如说，输入示例方程式，分别设 $C_{u}、HNO_{3}、Cu(NO_{3})2、NO、H_{2}O$ 的系数为 $x_{1}、x_{2}、x_{3}、x_{4}、x_{5}$ 。

由元素 $C_{u}、H、N、O$ 的守恒，得到方程组：

${x1=x3x2=2x5x2=2x3+x43x2=6x3+x4+x5\left\{\begin{matrix} x_{1} = x_{3}~~~~~~~~~~~~~~~~~~~~~ \\ x_{2} = 2x_{5}~~~~~~~~~~~~~~~~~~~ \\ x_{2} = 2x_{3}+x_{4} ~~~~~~~~~~~\\ 3x_{2} = 6x_{3}+x_{4}+x_{5} \end{matrix}\right.$

但方程组只有 $4$ 个方程，未知数却有 $5$ 个，无法求出唯一确定的一组解。

因为化学方程式系数间只是一种比例关系，可令 $x_{5}=1$ ，得：

${x1=34x2=2x3=34x4=12\left\{\begin{matrix} x_{1} = \frac{3}{4} \\ x_{2} = 2 \\ x_{3} = \frac{3}{4} \\ x_{4} = \frac{1}{2} \end{matrix}\right.$

各项系数都乘 $4$ ，得到最后的结果。

用计算机实现，应该用矩阵代替，消元法相同只是从操作方程变成了操作矩阵。

此外，还得将化学方程式转为线性方程组，难点在于去括号。

// 运行：在命令行输入 gcc 当前源文件.c #include <stdio.h> #include <string.h> #include <math.h> const double eps = 1e-10; int n, m, hash[ 1<<20 ]; double a[ 205 ] [ 205 ]; char z[ 205 ]; int anw, now;void Swap(double a, double b) {double t = a;a = b;b = t; } int get_number( int &i ) {int s = 0;while( '0' <= z[ i ] && '9' >= z[ i ] )s = ( s << 3 ) + ( s << 1 ) + z[ i ++ ] - '0';return s>1? s:1; }int get_str( int &i ) {if( z[ i ] > 'Z' || z[ i ] < 'A' )return -1;int s = z[ i ++ ];while( 'a' <= z[ i ] && z[ i ] <= 'z' )s = s * 10007 + z[ i ++ ];return s & ( ( 1<<20 )-1 ); }void countt( int l, int r, int f) {if( l == r ) return;int i = l, j = r;while( i < r - 1 && z[ i ] != '(' )i ++;while( j > l && z[ j ] != ')' )j --;if( z[ i ] == '(' ){countt(l, i, f);int w = j + 1, s = get_number(w);countt(i+1, j, f*s);countt(w, r, f);return;}for( i = l; i < r; ){int str = get_str(i);if( ! hash[ str ] )hash[ str ] = ++n;int hs = hash[ str ];a[ hs ] [ m ] += f * get_number(i);} }void init( ){printf("化学方程式配平前>> ");scanf("%s",z);int l = strlen(z), f = 1;z[ l ] = '#';for( int i = 0; i < l; ){int j = i;while( z[ j ] != '+' && z[ j ] != '=' && z[ j ] != '#' )j ++;m ++;countt(i, j, f);if( z[ j ] == '=' )f *= -1;i = j +1;} }void gs( ){ // 高斯消元for( int j = 1, i; j < m; j ++ ){for( i = now + 1; fabs(a[ i ] [ j ]) < eps && i <= n; i ++ );if( i > n ){anw = -1;continue;}now ++;for( int k = 1; k <= m; k ++ )Swap(a[i][k], a[now][k]);double ss = a[now][j];for( int k = 1; k <= m; k ++ )a[now][k] /= ss;for( i = 1; i <= n; i ++ )if( fabs(a[i][j]) > eps && i != now ){double s = a[i][j];for( int k = 1; k <= m; k ++ )a[i][k] -= a[now][k] * s;}} }int ans[205]; void solve( ){for( int i = now+1; i <= n; i ++ )if( fabs(a[i][m])>eps ){printf("无解\n");return;}if( anw == -1 ) puts("多组解法");for( int i = 1; i <= 1000; i ++ ){int p = 1;for( int j = 1; j < m; j ++ ){if( fabs(int(a[j][m] * i * -1 + 0.5) - a[j][m] * i * -1) > eps )p = 0;}if( p == 0 ) continue;for( int j = 1; j < m; j ++ )ans[j] = int( a[j][m] * i * -1 + 0.5 );ans[m] = i;printf("化学方程式配平后>> ");if( ans[1] != 1 ) printf("%d", ans[1]);int l = strlen(z), k = 1;for( int j = 0; j < l-1; j ++ ){putchar(z[j]);if( z[j] == '=' || z[j] == '+' )printf("%d",ans[++ k]);}break;} }int main(void){init( );gs( );solve( );return 0; }

逆矩阵求解

求解线性方程组除了消元法之外，逆矩阵也可求解。

在《矩阵实验：图形图像处理》，着重写了把矩阵看成一种对向量的函数、变换。

上面的方程，从线性系统的角度来看线性方程组的求解过程：

$[32−233−122−1]∗[x1x2x3]=[−4−54]\begin{bmatrix} 3 &2 &-2 \\ 3& 3 &-1 \\ 2& 2 &-1 \end{bmatrix}*\begin{bmatrix} x_{1}\\ x_{2}\\ x_{3} \end{bmatrix}=\begin{bmatrix} -4\\ -5\\ ~~4 \end{bmatrix}$

线性系统相当于：

$B * ? = y$

上面的方程，从线性变换的角度来看线性方程组的求解过程：

$[32−233−122−1]∗[−3021−22]=[y1y2y3]\begin{bmatrix} 3 &2 &-2 \\ 3& 3 &-1 \\ 2& 2 &-1 \end{bmatrix}* \begin{bmatrix} -30\\ ~~21\\ -22\end{bmatrix} =\begin{bmatrix} y_{1}\\ y_{2}\\ y_{3} \end{bmatrix}$

线性变换相当于；

$B * x = ?$ ，矩阵 $B$ 是变换。

看线性系统，已知变换矩阵 $B$ 、输出矩阵 $y$ ，求解输入矩阵 $x$ 。

那求解线性方程组就是一个逆变换的过程；因此，有一个逆矩阵表征这种变换。

逆矩阵：若 $A B = E$ ，则称 $A, B$ 互为逆矩阵；记作： $B^{-1}=A ~~or~~A^{-1}=B$ 。

如果我们在等式俩边都乘，式子如下：

$Ex=(AB)x=A(Bx)=B^{-1}(Bx)$

这个式子在几何上，变换来变换去，最后回到初始状态：

参照这种思想(矩阵是一种变换的思想)， $B * ? = y$ 可以用 $x=B^{-1}*y$ 来求！！

求线性方程组就变成了求解某个矩阵的逆矩阵。

现在，唯一要明白的是：逆矩阵怎么求？

矩阵的逆

除零以外，一个数字乘以这个数字的倒数（逆）等于一：

$x*x^{-1} = 1$

这是在数字系统里，事实上，矩阵也可以逆：

$A B = B A = I$ ， $I$ 是单位矩阵，相当于数字里面的 $1$ 。

如果矩阵满足上述 $等式$ ，则称 $B$ 是 $A$ 的逆矩阵，记做： $B=A^{-1}$

大部分矩阵都是有逆的：

有逆的矩阵叫可逆矩阵 or 非奇异矩阵
不可逆的矩阵叫不可逆矩阵 or 奇异矩阵

因为矩阵乘法不满足交换律，所以可逆分成了俩种：

左逆： $B A = I$ ，只有左逆，是不可逆矩阵；
右逆： $A B = I$ ，只有右逆，是不可逆矩阵；

可逆矩阵是同时有左逆和右逆，只有左逆或只有右逆都不叫可逆。

只有方阵可逆（行数 = 列数），因为对于一个方阵来说，左逆和右逆一定是同时存在的。

本科的线性代数主要研究方阵（除了线性系统）。

矩阵求逆：

# 运行：在命令行里输入 python 当前源文件.py import numpy as np from scipy import linalgA = np.array([[1, 35, 0],[0, 2, 3],[0, 0, 4]])A_n = linalg.inv(A) print(A_n) print(np.dot(A, A_n))

运行结果：

[[ 1. -17.5 13.125][ 0. 0.5 -0.375][ 0. 0. 0.25 ]][[ 1. 0. 0.][ 0. 1. 0.][ 0. 0. 1.]]

C++ 也有专属的线代库 armadillo，因为有数据类型有精度限制可以用 Boost.Multiprecision，用于科学计算也就很方便。

解的结构

消元法所做的是把一个增广矩阵变换为一个阶梯型矩阵。

阶梯型矩阵：

非零行第一个元素（主元）为 $1$
主元所在列的其他元素均为 $0$

满足这种条件的矩阵，也被称为 “行最简形式”。

线性方程组通过消元法变成行最简形式，这样就可以直接判断方程组解的情况：

适定方程组：方程组有唯一解，理想情况，遇到的情况不多；
超定方程组：方程组无解，无论三个未知数取什么解，方程始终不能满足，因此方程组无解；遇到的概率很大；
欠定方程组：方程组有无数解，这种方程组的未知数有无数个解，因此实用价值不大。

一般的判定方法，是通过矩阵的秩来判断方程组的类型。

欠定方程组（无数解）：系数矩阵的秩 = 增广矩阵的秩 < 未知数个数；
适定方程组（唯一解）：系数矩阵的秩 = 增广矩阵的秩 = 未知数个数；
超定方程组（无解）：系数矩阵的秩 $≠\neq$ 增广矩阵的秩。

超定方程组近似估值

我们遇到的多是 超定方程组，因此学习一下如何求超定方程的近似解。

先看一个比较简单的问题。

设 $→u=(1,1)T,→v=(2,0)T,→b=(3,1)T\underset{u}{\rightarrow}=(1,1)^{T},\underset{v}{\rightarrow}=(2,0)^{T},\underset{b}{\rightarrow}=(3,1)^{T}$ ，能否找到一组 $x_{1}, x_{2}$ 满足： $x1∗→u+x2∗→v=→bx_{1}*\underset{u}{\rightarrow}+x_{2}*\underset{v}{\rightarrow}=\underset{b}{\rightarrow}$ 。

用矩阵表示： $[u,v]∗[x1x2]=b\left [ u,v \right ]*\begin{bmatrix} x_{1}\\ x_{2} \end{bmatrix}=b$ ，也等同于 $[1210]∗[x1x2]=[31]\begin{bmatrix} 1 &2 \\ 1 & 0 \end{bmatrix}*\begin{bmatrix} x_{1}\\ x_{2} \end{bmatrix}=\begin{bmatrix} 3\\ 1 \end{bmatrix}$ ， (解出来是，适定方程组)，解得(唯一解) ： $x_{1}=1, x_{2}=1$ 。

第二个问题，设 $→u=(1,1,0)T,→v=(2,0,0)T,→b=(3,1,3)T\underset{u}{\rightarrow}=(1,1,0)^{T},\underset{v}{\rightarrow}=(2,0,0)^{T},\underset{b}{\rightarrow}=(3,1,3)^{T}$ ，试问，能否找到一组 $x_{1}, x_{2}$ 满足： $x1∗→u+x2∗→v=→bx_{1}*\underset{u}{\rightarrow}+x_{2}*\underset{v}{\rightarrow}=\underset{b}{\rightarrow}$ 。

用矩阵表示： $[u,v]∗[x1x2]=b\left [ u,v \right ]*\begin{bmatrix} x_{1}\\ x_{2} \end{bmatrix}=b$ ，也等同于 $[121000]∗[x1x2]=[313]\begin{bmatrix} 1 &2 \\ 1 &0 \\ 0 & 0 \end{bmatrix}*\begin{bmatrix} x_{1}\\ x_{2} \end{bmatrix}=\begin{bmatrix} 3\\ 1\\ 3 \end{bmatrix}$ ，(解出来是，超定方程组)，无解。

图中的向量 $b$ ，脱离了 $u, v$ 俩个向量的平面；这就从几何上说满足条件的未知数是无法找到的。

虽然无解，但却可以找到一个最接近向量 $b$ 的 $b_{1}$ 向量，如下图。

$b$ 向量 - $b_{1}$ 向量 = $e$ 向量

几何推导过程：

方程组表达式也可以推导：

方程组表达式： $A→x=→bA\underset{x}{\rightarrow}=\underset{b}{\rightarrow}$
在等式俩边都乘一个 $AT：ATA→x=AT→bA^{T}：A^{T}A\underset{x}{\rightarrow}=A^{T}\underset{b}{\rightarrow}$

$A$ 的转置乘 $A$ 必然是一个实对称矩阵，实对称矩阵必定可逆，等式的左右俩边还可以乘 $A$ 的转置： $(ATA)−1AT→b(A^{T}A)^{-1}A^{T}\underset{b}{\rightarrow}$ 。

工程应用：天体轨道参数估计

前置知识：《极坐标、开普勒第一定律》。

根据开普勒第一定律，当忽略其他天体的重量吸引时，一个天体应该取椭圆、抛物线或双曲线轨道。

在极坐标 $(Θ,r)(\Theta,~r)$ 中，天体的位置满足一个方程：

$\beta + e(r~cos(\Theta))$

$β\beta$ 为轨道常数， $e$ 是轨道偏心率；对于椭圆 $0 < e < 1$ ，对于抛物线 $e > 1$ 。

著名的哈雷彗星是一个椭圆，台湾的鹿林彗星是一个抛物线。

天文学家利用太空望远镜🔭，观测到一个新的天体，上图的原点 $O$ 代表地球，坐标是极坐标系。

新天体，第一次观测出现的位置（得到极坐标 $(Θ,r)(\Theta,~r)$ ：）：

连续观测了五次，得到了五组数据（五个坐标）：

而后，要做的就是根据这些已有的数据推导天体的轨道方程！！

也就是，估计轨道常数参数 $β\beta$ 、轨道偏心率 $e$ 。

因为我们已经得到了五组 $(Θ,r)(\Theta,~r)$ ，而天体的位置又必然满足轨道方程 $\beta + e(r~cos(\Theta))$ ，就形成了五组二元一次方程。

${1⋅β+3.00cos(0.88)⋅e=3.001⋅β+2.30cos(1.10)⋅e=2.301⋅β+1.65cos(1.42)⋅e=1.651⋅β+1.77cos(1.25)⋅e=1.251⋅β+2.14cos(1.01)⋅e=1.01\left\{\begin{matrix} 1· \beta + 3.00~cos(0.88) ·e = 3.00 \\ 1· \beta + 2.30~cos(1.10) ·e = 2.30 \\ 1· \beta + 1.65~cos(1.42) ·e = 1.65 \\ 1· \beta + 1.77~cos(1.25) ·e = 1.25 \\ 1· \beta + 2.14~cos(1.01) ·e = 1.01 \end{matrix}\right.$

看这些方程的类型，就是超定方程组 — 系数还有小数，基本消不了。

用矩阵表示（矩阵乘法）：

$[13.00cos(0.88)12.30cos(1.10)11.65cos(1.42)11.77cos(1.25)12.14cos(1.01)]∗[βe]=[3.002.301.651.251.01]\left[ \begin{matrix} 1 & 3.00~cos(0.88)\\ 1 & 2.30~cos(1.10)\\ 1 & 1.65~cos(1.42)\\ 1 & 1.77~cos(1.25)\\ 1 & 2.14~cos(1.01) \end{matrix} \right]* \left[ \begin{matrix} \beta\\ e \end{matrix} \right]= \left[ \begin{matrix} 3.00\\ 2.30\\ 1.65\\ 1.25\\ 1.01 \end{matrix} \right]$

解得： $β=1.45,e=0.81\beta = 1.45, ~e = 0.81$ 。

所以，天体轨道方程为： $1.45+0.81r~cos(\Theta)$ 。

整理得： $r=1.451−0.81cos(Θ)r=\frac{1.45}{1-0.81~cos(\Theta)}$

观测的数据越多，参数估计就越准确。

# 运行：在命令行输入 python 当前源文件.py import numpy as np import matplotlib.pyplot as plt# 1.准备好观测到的 (θ,r) 数据 r = 贝塔 + e(r cosθ) data = np.array([ [ 0.88, 3.00],[ 1.10, 2.30],[ 1.42, 1.65],[ 1.77, 1.25],[ 2.14, 1.01] ])# 2.构建系数矩阵A A = np.ones((data.shape[0],data.shape[1])) for i in range(0,data.shape[0]): # i从 [0,5) 区间逐一取值theta , r = data[i][0] , data[i][1]A[i][1] = r * np.cos(theta) print("Mat A:\n", A)# 3.构造常数项矩阵b b = data[:,1] print("Mat b:\n", b)# 4. 根据估值公式求解参数 x = np.dot(np.dot(np.linalg.inv(np.dot(A.T,A)), A.T),b) print("--------------\nMat x:\n",x)# 5. 画出天体运行轨道 theta_ = np.linspace(0, 2*np.pi, 100) r_ = x[0] / (1-x[1]* np.cos(theta_)) graph = plt.subplot(111, polar=True) # 1行1列图像中的第一个 data_= data.T graph.plot(data_[0,:],data_[1,:],'go') # 将已测量数据打上绿色的o（'go'） graph.plot(theta_, r_ ,'r', linewidth=1) plt.show()

效果演示：

总结

以上是生活随笔为你收集整理的线性系统实验：化学方程式配平与天体轨道参数估计的全部内容，希望文章能够帮你解决所遇到的问题。

如果觉得生活随笔网站内容还不错，欢迎将生活随笔推荐给好友。

上一篇：里恩临床试验电子数据采集系统(EDC)介
下一篇：光盘在服务器上装系统,如何通过DELL引