求解方程( 迭代法 牛顿迭代法 二分法)

时间:2022-10-27 18:29:37
迭代的意思是反反复复地执行某一步骤、程序或者事件,在数学中用的比较常见。

以上代码转自https://blog.csdn.net/pengwill97/article/details/77200372

https://blog.csdn.net/akatsuki__itachi/article/details/80719686

首先,迭代法解方程的实质是按照下列步骤构造一个序列x0,x1,…,xn,来逐步逼近方程f(x)=0的解:

1)选取适当的初值x0;

2)确定迭代格式,即建立迭代关系,需要将方程f(x)=0改写为x=φ(x)的等价形式;

3)   构造序列x0,x1,……,xn,即先求得x1=φ(x0),再求x2=φ(x1),……如此反复迭代,就得到一个数列x0, x1,……,xn,若这个数列收敛,即存在极值,且函数 φ(x)连续,则很容易得到这个极限值求解方程( 迭代法 牛顿迭代法 二分法),x*就是方程f(x)=0的根。

#define eps 1e-8
int main()
{
    x0=初始近似根;
    do{
        x1=x0;
        x0=g(x1); //按特定的方程计算新的近似根
    }while(fabs(x0-x1)>Epsilon);
    printf("方程的近似根是%f\n",x0);
}

 注意:如果方程无解,算法求出的近似根序列就不会收敛,那么迭代过程就会变成死循环。因此,在使用迭代算法前应先考察方程是否有解,并在算法中对迭代次数给予限制。

 

#include<stdio.h>
#include<iostream>
#include<math.h>
#define eps 1e-8
using namespace std;
void NewtonIteration();//牛顿迭代法求解方程
void DichotomySolving();//二分法求解方程
int main ()
{
    NewtonIteration();
    DichotomySolving();
}
//牛顿迭代法:xn=x(n-1)-f(x(n-1))/f'(x(n-1))
void NewtonIteration()
{
    int a = 3, b = 2, c = 1, d = -6;//系数
    float x1 = 1, x0;
    float f0, f1;//f0是原方程 f1是方程的一阶导
    do
    {
        x0 = x1;
        f0 = ((a * x0 + b) * x0 + c) * x0 + d;//即a*x^3+b*x^2+c*x+d
        f1 = (3 * a * x0 + 2 * b) * x0 + c;//一阶导数
        x1 = x0 - f0/f1;
    }while(fabs(x1 - x0) >= eps);
    printf("(%d)x^3+(%d)x^2+(%d)x+(%d) = 0\n", a, b, c, d);
    cout << "方程的解为:" << x1 << endl;
}
 
void DichotomySolving()//二分法求解方程
{
    float x, x1 = 0, x2 = 2;
    float f1, f2, f;
    cout << "input x1, x2 (f(x1)*f(x2)<0)" << endl;
    //假设原方程为 1/2x^3+2x^2-8
    f1 = pow(x1, 3) / 2 + pow(x1, 2) * 2 - 8;
    f2 = pow(x2, 3) / 2 + pow(x2, 2) * 2 - 8;
    if(f1 * f2 > 0)
    {
        cout << "error" << endl;
        return;
    }
    do
    {
        x = (x1 + x2) / 2;
        f = pow(x, 3) / 2 + 2 * pow(x, 2) - 8;//f为飞f1
        if(f == 0)
            break;
        if(f1 * f > 0)//已知解一定在f1与f2之间 因此以f1为切入点逐步逼近解
        {//f1 * f > 0证明解在右半部分中 因此要

高斯消元法

#include<stdio.h>
#include<algorithm>
#include<iostream>
#include<string.h>
#include<math.h>
using namespace std;
const int MAXN=50;
int a[MAXN][MAXN];//增广矩阵
int x[MAXN];//解集
bool free_x[MAXN];//标记是否是不确定的变元
int gcd(int a,int b)
{
    if(b == 0) return a;
    else return gcd(b,a%b);
}
inline int lcm(int a,int b)
{
    return a/gcd(a,b)*b;//先除后乘防溢出
}
// 高斯消元法解方程组(Gauss-Jordan elimination).(-2表示有浮点数解,但无整数解,
//-1表示无解,0表示唯一解,大于0表示无穷解,并返回*变元的个数)
//有equ个方程,var个变元。增广矩阵行数为equ,分别为0到equ-1,列数为var+1,分别为0到var.
int Gauss(int equ,int var)
{
    int i,j,k;
    int max_r;// 当前这列绝对值最大的行.
    int col;//当前处理的列
    int ta,tb;
    int LCM;
    int temp;
    int free_x_num;
    int free_index;

    for(int i=0; i<=var; i++)
    {
        x[i]=0;
        free_x[i]=true;
    }
    //转换为阶梯阵.
    col=0; // 当前处理的列
    for(k = 0; k < equ && col < var; k++,col++) // 枚举当前处理的行.
    {
        // 找到该col列元素绝对值最大的那行与第k行交换.(为了在除法时减小误差)
        max_r=k;
        for(i=k+1; i<equ; i++)
        {
            if(abs(a[i][col])>abs(a[max_r][col])) max_r=i;
        }
        if(max_r!=k) // 与第k行交换.
        {
            for(j=k; j<var+1; j++) swap(a[k][j],a[max_r][j]);
        }
        if(a[k][col]==0) // 说明该col列第k行以下全是0了,则处理当前行的下一列.
        {
            k--;
            continue;
        }
        for(i=k+1; i<equ; i++) // 枚举要删去的行.
        {
            if(a[i][col]!=0)
            {
                LCM = lcm(abs(a[i][col]),abs(a[k][col]));
                ta = LCM/abs(a[i][col]);
                tb = LCM/abs(a[k][col]);
                if(a[i][col]*a[k][col]<0)tb=-tb;//异号的情况是相加
                for(j=col; j<var+1; j++)
                {
                    a[i][j] = a[i][j]*ta-a[k][j]*tb;
                }
            }
        }
    }
    // 1. 无解的情况: 化简的增广阵中存在(0, 0, ..., a)这样的行(a != 0).
    for (i = k; i < equ; i++)  // 对于无穷解来说,如果要判断哪些是*变元,那么初等行变换中的交换就会影响,则要记录交换.
    {
        if (a[i][col] != 0) return -1;
    }
    // 2. 无穷解的情况: 在var * (var + 1)的增广阵中出现(0, 0, ..., 0)这样的行,即说明没有形成严格的上三角阵.
    // 且出现的行数即为*变元的个数.
    if (k < var)
    {
        return var - k; // *变元有var - k个.
    }
    // 3. 唯一解的情况: 在var * (var + 1)的增广阵中形成严格的上三角阵.
    // 计算出Xn-1, Xn-2 ... X0.
    for (i = var - 1; i >= 0; i--)
    {
        temp = a[i][var];
        for (j = i + 1; j < var; j++)
        {
            if (a[i][j] != 0) temp -= a[i][j] * x[j];
        }
        if (temp % a[i][i] != 0) return -2; // 说明有浮点数解,但无整数解.
        x[i] = temp / a[i][i];
    }
    return 0;
}
int main(void)
{
    //    freopen("in.txt", "r", stdin);
    //    freopen("out.txt","w",stdout);
    int i, j;
    int equ,var;
    while (scanf("%d %d", &equ, &var) != EOF)
    {
        memset(a, 0, sizeof(a));
        for (i = 0; i < equ; i++)
        {