另一个缘由是封锁不够,问这串数字符合f

题意:给你一串数字,问那串数字符合f[n]
= a*f[n-1],f[n] = a*f[n-1]+b*f[n-2],f[n] =
a*f[n-1]+b*f[n-2]+c*f[n-3]这个方程中的哪个,然后要你提交第n+1项,假若符合四个方程,项数小的先行(第一个方程优先)。

下边是运行一个adams/car模型出现的谬误。

高斯消元的真面目就是模拟解方程

解法:这题我先处理看是否满意f[n] =
a*f[n-1]的花样,如若不满足,则用高斯消元借出两项和三项的境况的a,b,c,比如第二个方程,f[3]
= a*f[2]+b*f[1],f[4] =
a*f[3]+b*f[2],三个方程两个未知量,用高斯消元解出a,b,这里恐怕不是整数,我将他们加了个0.5取下整,居然对了。后来看这场交锋没一个人是用的高斯消元,所以不领会这么是不是科学,有看出来端倪的迎接评论告诉自己。

—- ERROR —-
   The system matrix has a zero pivot for column 2142, which is
associated
   with VARIABLE/64 Algb Var.  Consequently, the matrix is numerically
singular.

想象一下,你平日解n元三遍方程组的时候是咋办的?

代码:

谬误分析

答案是逐日消元啦~

 

我们可能都知情当出现zero
pivot错误的时候,一般是model中冒出了过约束(over
constraint)或者是束缚不够而致使rigid mode.
以前自己是把这看作一条理论来记得,原因相比较模糊,前天把它搞搞精通,落在文字上与我们大快朵颐。说pivot,
得先从高斯消元讲起。有限元软件求解刚度矩阵时,一般是用高斯消元。对于高斯消元,我们应该相比较熟谙,这是一种基于先正向消元,再反向迭代求解的求解办法(It
is based on triangularization of the coefficient matrix and evaluation
of the unknowns by back-subsitution starting from the last
equation)。高斯消元后的矩阵是这般的

对此方程组:

710官方网站 1710官方网站 2

| a11  a12   a13  a14… a1n  |      x1       c1

a11*x1+a12*x2+a13*x3+……+a1n*xn=b1

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
#define N 4

int f[14];
typedef double Matrix[N][N];
int x,y,z;

void gauss_elimination(Matrix A,int n)
{
    int i,j,k,r;
    for(i=0;i<n;i++)
    {
        //选一行r并与i行交换
        r = i;
        for(j=i+1;j<n;j++)
            if(fabs(A[j][i]) > fabs(A[r][i]))
                r = j;
        if(r != i)
        {
            for(j=0;j<=n;j++)
                swap(A[r][j],A[i][j]);
        }
        //与第i+1~n行进行消元
        for(k=i+1;k<n;k++)
        {
            double f = A[k][i]/A[i][i];  //为了让A[k][i] = 0,第i行乘以的倍数
            for(j=i;j<=n;j++)
                A[k][j] -= f*A[i][j];
        }
    }
    //回代
    for(i=n-1;i>=0;i--)
    {
        for(j=i+1;j<n;j++)
            A[i][n] -= A[j][n]*A[i][j];
        A[i][n] /= A[i][i];
    }
    x = (int)floor(A[0][n]+0.5);
    y = (int)floor(A[1][n]+0.5);
    if(n == 3)
        z = (int)floor(A[2][n]+0.5);
}

int main()
{
    int t,n,i,j;
    scanf("%d",&t);
    while(t--)
    {
        scanf("%d",&n);
        memset(f,0,sizeof(f));
        for(i=1;i<=n;i++)
            scanf("%d",&f[i]);
        int ans = Mod;
        int a1,a2,a3;
        int flag;
        if((f[1] == 0 && f[2] == 0) || f[2]%f[1] == 0)
        {
            if(f[1] == 0 && f[2] == 0)
                a1 = 1;
            else
                a1 = f[2]/f[1];
            flag = 1;
            for(i=3;i<=n;i++)
            {
                if(f[i] != a1*f[i-1])
                    flag = 0;
            }
            if(flag)
                ans = a1*f[n];
        }
        if(ans != Mod)
        {
            printf("%d\n",ans);
            continue;
        }
        Matrix A;
        A[0][0] = A[1][1] = f[2];
        A[0][1] = f[1];
        A[1][0] = f[3];
        A[0][2] = f[3];
        A[1][2] = f[4];
        gauss_elimination(A,2);
        flag = 1;
        for(i=3;i<=n;i++)
        {
            if(f[i] != x*f[i-1]+y*f[i-2])
                flag = 0;
        }
        if(flag)
            ans = x*f[n]+y*f[n-1];
        if(ans != Mod)
        {
            printf("%d\n",ans);
            continue;
        }
        A[0][0] = A[1][1] = A[2][2] = f[3];
        A[0][1] = A[1][2] = f[2];
        A[0][2] = f[1];
        A[1][0] = A[2][1] = f[4];
        A[2][0] = f[5];
        A[0][3] = f[4];
        A[1][3] = f[5];
        A[2][3] = f[6];
        gauss_elimination(A,3);
        //printf("%d %d %d\n",x,y,z);
        ans = x*f[n]+y*f[n-1]+z*f[n-2];
        if(ans != Mod)
            printf("%d\n",ans);
    }
    return 0;
}

| 0      a22   a23  a24… a2n  |      x2       c2

a21*x1+a22*x2+a23*x3+……+a2n*xn=b2

View Code

| 0      0       a33  a34…a3n   |      x3       c3

a31*x1+a32*x2+a33*x3+……+a3n*xn=b3

 

| 0      0       0      a44…a4n  |   {  x4}={c4}

……..

| ………………………… …….|      …       …

formula1:把x1周全化为1,即方程除以a11

| 0      0        0        0 … ann |       xn      cn

然后对于剩下的n-1个方程组,用剩下的n-1个方程组去减formula1*ai1【ai1为相应的方程组的首先个全面】

诸如此类第一步先求出xn,第二步就足以求出x(n-1),如此类推。在这一个矩阵中每行的首先个非零的系数就是pivot。那么zero
pivot的意思就明确了,就是指在高斯消元后的刚度矩阵中现身了一个全为零的一条龙。一个原因是出新了过构束,就好比去用10个方程去解一个9个未知数,一定有一个方程可以消去。多余这几个方程有可能与原来的某一方程等价,也有可能与某一方程冲突,但结果都是zero
pivot。 另一个缘由是约束不够,有力,但与之对应没有刚度,无疑会出现zero
pivot。

这规范剩余的每个方程的x1就被消去啊

一旦是束缚不够时,一般在message file里还会出现NUMERICAL SINGULARITY
warnings。这相似是因为力除以0刚度出现了无穷大的活动。

下一场用formula2去消掉x2,渐渐消掉每个x,直至最后一个方程只剩an,算出an,往回迭代,就解出这么些方程了~

如上就是高斯消元的有血有肉考虑步骤。

要留意的是,在除一个数从前若发现那么些数为0,那么无解

来道水题练练手= =

1013: [JSOI2008]球形空间发出器sphere

710官方网站,Time Limit: 1 Sec  Memory
Limit: 162 MB
Submit: 5976  Solved: 3115
[Submit][Status][Discuss]

Description

  有一个球形空间暴发器能够在n维上空中暴发一个僵硬的球体。现在,你被困在了这一个n维球体中,你只知道球
表面n+1个点的坐标,你需要以最快的进度规定这些n维球体的球心坐标,以便于摧毁这多少个球形空间发生器。

Input

  第一行是一个整数n(1<=N=10)。接下来的n+1行,每行有n个实数,表示球面上一点的n维坐标。每一个实数精确到小数点
后6位,且其相对值都不超越20000。

Output

  有且只有一行,依次给出球心的n维坐标(n个实数),两个实数之间用一个空格隔开。每个实数精确到小数点
后3位。数据保证有解。你的答案必须和业内输出一模一样才可以得分。

Sample Input

2
0.0 0.0
-1.0 1.0
1.0 0.0

Sample Output

0.500
1.500

HINT

  提醒:给出六个概念:1、 球心:到球面上无限制一点离开都分外的点。2、
距离:设五个n为空间上的点A, B

的坐标为(a1, a2, …, an), (b1, b2, …, bn),则AB的距离定义为:dist = sqrt(
(a1-b1)^2 + (a2-b2)^2 + 

… + (an-bn)^2 )

俺们了然,圆心到圆上所有点的离开相等,这样子我们就足以列出n个方程,两边平方去括号可以发现xi^2都可以约去,最后只剩一个这么的方程组:

2*(a(i)1-a(i-1)1)x1+2*(a(i)x-a(i-1)x)xx+……+2*(a(i)n-a(i-1)1n)xn=a(i)1^1+a(i)2^2+……-a(i-1)1^2-a(i-1)2^2-……..

简直就是板题呐

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define LL long long int
using namespace std;
const int maxn=15,INF=2000000000,P=1000000007;

double A[maxn][maxn],num[maxn][maxn],ans[maxn];
int n;

void init(){
    cin>>n;
    for(int i=0;i<=n;i++)
        for(int j=1;j<=n;j++)
            scanf("%lf",&num[i][j]);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++){
            A[i][j]=2*num[i][j]-2*num[i-1][j];
            A[i][n+1]+=num[i][j]*num[i][j]-num[i-1][j]*num[i-1][j];
        }
}

void gaosi(){
    for(int i=1;i<=n;i++){
        double t=A[i][i];
        if(t==0.0) return;
        for(int j=i;j<=n+1;j++)
            A[i][j]/=t;
        for(int j=i+1;j<=n;j++){
            t=A[j][i];
            for(int k=i;k<=n+1;k++)
                A[j][k]-=A[i][k]*t;
        }
    }
    for(int i=n;i>0;i--){
        for(int j=i+1;j<=n;j++)
            A[i][n+1]-=ans[j]*A[i][j];
        ans[i]=A[i][n+1]/A[i][i];
    }
}

void print(){
    printf("%.3lf",ans[1]);
    for(int i=2;i<=n;i++)
        printf(" %.3lf",ans[i]);
}

int main(){
    init();
    gaosi();
    print();
    return 0;
}

相关文章