常用计算方法(C语言代码)(计算方法课程)

目录

  • 矩阵相乘
  • 各种求解一元线性方程解的解法
  • 普通高斯消元
  • 列主元素高斯消元
  • 普通平方根法求解对称正定矩阵
  • 改进的平方根求法解对称正定矩阵
  • 矩阵三角分解(LU)(杜立特分解)
  • 追赶法 求解 三对角方程组
  • 雅可比迭代法
  • 高斯赛德尔迭代法及超松弛迭代法
  • 拉格朗日插值
  • 最小二乘法求解拟合曲线
  • 牛顿插值多项式 (均差)
  • 变步长梯形求积法
  • 复化梯形公式 与 复化辛普森公式求积
  • 改进的欧拉法 求解 微分方程的 数值
  • 数值积分 之 龙贝格求积分法
  • 常微分方程初值问题 数值解法 之 经典 - 四阶龙格库塔 法

说明

这些计算方法是我在计算方法课程实验时侯写的,之前放在github(喜欢的可以给star(heihei))代码仓库,现在搬到博客,方便以后自己复习,由于很多过程设计到公式,并且比较繁琐,就没有贴上详细的算法解释,到时候不理解可以看一下计算方法的书(我的是刘师少版)回顾一下。

矩阵相乘

这个比较简单,只需要对应行和列相称即可,一个三重循环。

//计算方法实验:矩阵相乘以及文件输入输出

#include <iostream>
#include <stdio.h>
#include <stdlib.h>

using namespace std;
const int M = 2;
const int N = 2;
const int L = 2;

int a[M][L],b[L][N],c[M][N];

int main(){
    FILE *in,*out;
    if((in = fopen("in.txt","r")) == NULL){cout<<"Cannot open the file"<<endl;exit(0);}
    if((out = fopen("out.txt","w")) == NULL){cout<<"Cannot open the file"<<endl;exit(0);}
    for(int i = 0; i < M; i++){
        for(int j = 0; j < L; j++){
            //cin>>a[i][j];
            fscanf(in,"%d",&a[i][j]);
        }
    }
    for(int i = 0; i < L; i++){
        for(int j = 0; j < N; j++){
            //cin>>b[i][j];
            fscanf(in,"%d",&b[i][j]);
        }
    }
    for(int i = 0; i < M; i++){
        for(int j = 0; j < N; j++){
            c[i][j] = 0;
            for(int k = 0; k < L; k++)c[i][j] += a[i][k] * b[k][j];
            cout<<c[i][j]<<" ";
            fprintf(out,"%d ",c[i][j]);
        }
        cout<<endl;
        fprintf(out,"\n");
    }
    return 0;
}

in.txt,out.txt文件示例

in.txt    out.txt
1 2       19 22           
3 4       43 50 
5 6
7 8

各种求解一元线性方程解的解法

一元线性方程求解的方法有很多,包括二分法收敛,简单迭代法,史蒂芬森加速,牛顿迭代法,牛顿下山法,弦截法,各种算法的原理不一一解说,直接上之前写过的代码:

//各种求解一元线性方程的数值解法总结:

#include <iostream>
#include <stdio.h>
#include <math.h>

using namespace std;
//f(x) = x*x*x - 3*x - 1有三个实根,x1 = 1.8793 ,x2 = -0.34727,x3 = -1.53209
//下面采用六种不同的计算格式,求解f(x) = 0的根 x1 或 x2
double f1(double x) { return (3 * x + 1) / (x*x); }
double fd1(double x) { return (3 * x*x + 6 * x*x + 2 * x) / (x*x*x*x); }//f1的导数
double f2(double x) { return (x*x*x - 1) / 3; }
double fd2(double x) { return x*x; }//f2的导数
double f3(double x) { return pow(3 * x + 1, 1.0 / 3); }
double f4(double x) { return 1 / (x*x - 3); }
double fd4(double x) { return 2 * x / ((x*x - 2)*(x*x - 2)); }
double f5(double x) { return sqrt(3 + 1 / x); }
double f6(double x) { return x - (x*x*x - 3 * x - 1) / (3 * (x*x - 1)); }
//f7是书上弦截法的例题
double f7(double x) { return exp(-x) - x; }
double fd7(double x) { return -exp(-x) - 1; }
double eps = 1e-7;

int erFen(double x0, double x1) {
    int sum = 0;
    if (f1(x0) * f1(x1) < 0) {
        do {
            double mid = (x0 + x1) / 2;
            if (f1(mid)*f1(x0) < 0)x1 = mid;
            else x0 = mid;
            sum++;
            printf("x = %8.6f\n", x1);
        } while (fabs(x1 - x0) > eps);
    }
    return sum;
}

int Simple(double x0, double x1) {//简单迭代法
    double d = 0;
    int sum = 0;
    do {
        x1 = f2(x0);
        d = fabs(x1 - x0);
        x0 = x1;
        sum++;
        printf("x = %8.6f\n", x1);
    } while (d >= eps);
    return sum;
}

int SiTiFen(double x0, double x1) {//史蒂芬森加速算法
    double d = 0, y0 = 0, z0 = 0;
    int sum = 0;
    do{
        y0 = f4(x0);
        z0 = f4(y0);
        x1 = z0 - (z0 - y0)*(z0 - y0) / (z0 - 2 * y0 + x0);
        d = fabs(x1 - x0);
        x0 = x1;
        sum++;
        printf("x = %8.6f\n", x1);
    } while (d >= eps);
    return sum;
}

int NiuDun(double x0, double x1) {//牛顿跌代法
    double d = 0;
    int sum = 0;
    do {
        x1 = x0 - f2(x0) / fd2(x0);
        d = fabs(x1 - x0);
        x0 = x1;
        sum++;
        printf("x = %8.6f\n", x1);
    } while (d >= eps);
    return sum;
}

int XiaShan(double x0,double x1){//牛顿下山法
    double d = 0;
    int sum = 0;
    do {
        if (fd2(x0) == 0) break; //倒数为0
        double t = 1.0;
        while (1) {
            x1 = x0 - (1.0 / t)*(f2(x0) / fd2(x0));
            if (fabs(f2(x1)) < fabs(f2(x0)))break;
            else {
                if (t == 1.0)t = t + 1;
                else t = t + 2;  //取1,1/2,1/4,1/8
            }
        }
        d = fabs(x1 - x0);
        sum++;
        x0 = x1;
        printf("x = %8.6f\n", x1);
    } while (d >= eps);
    return sum;
}

int XianJie(double x0, double x1, double x2) {//弦截法
    //答案:解:0.56714  (f7)
    int sum = 0;
    do {
        x2 = x1 - f7(x1) * (x1 - x0) / (f7(x1) - f7(x0));
        x0 = x1;
        x1 = x2;
        sum++;
        printf("x = %8.6f\n", x2);
    } while (fabs(x1 - x0) >= eps);
    return sum;
}

int main() {
    double x0 = 0,x1 = 0, x2 = 0;
    x0 = -0.5, x1 = -0.25;//选取二分的开始两个点
    printf("二分法收敛次数是 %d\n\n", erFen(x0, x1));
    x1 = 0, x0 = 0.5;
    printf("简单迭代法收敛次数是 %d\n\n", Simple(x0, x1));
    x1 = 0, x0 = 0.5;
    printf("史蒂芬森加速算法收敛次数是 %d\n\n", SiTiFen(x0, x1));
    x1 = 0, x0 = 0.5;
    printf("牛顿迭代法收敛次数是 %d\n\n", NiuDun(x0, x1));//牛顿法和下山法貌似没有算出正确答案
    x1 = 0, x0 = 0.5;
    printf("牛顿下山法收敛次数是 %d\n\n", XiaShan(x0, x1));
    x1 = 0.5, x0 = 0.6;
    printf("弦截法收敛次数是 %d\n\n", XianJie(x0, x1, x2));
    return 0;
}

普通高斯消元

高斯消元求解线性方程组的解是最经典的解法。也不多说,直接上代码,其中包括每一次消元之后的增广矩阵

#include<iostream>
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <cmath>
#define N 1000
using namespace std;
double a[N][N],b[N],x[N];
int main(){
    FILE *in;
    int n;
    if((in = fopen("in.txt","r")) == NULL){ cout<<"cannot open the file"<<endl; exit(0);}
    fscanf(in,"%d",&n);
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n; j++)
            fscanf(in,"%lf",&a[i][j]);
        fscanf(in,"%lf",&b[i]);
    }
    int k;
    for(k = 1; k < n; k++){
        for(int i = k+1; i <= n; i++){
            for(int j = k+1; j <= n; j++)  ///注意一定要从 K + 1开始
                a[i][j] = a[i][j]-a[k][j]*a[i][k]/a[k][k];
            b[i] = b[i]-b[k]*a[i][k]/a[k][k];
        }
//        printf("经过第%d步后增广矩阵为:\n",k);
//        for(int l = 1; l <= n; l++){
//            for(int j = 1; j <= n; j++) printf("%f ",a[l][j]);
//            printf("%f ",b[l]);
//            printf("\n");
//        }
    }
    for(int i = n; i > 0; i--){
        double sum = 0;
        for(int j = i+1; j <= n; j++) sum = sum+a[i][j]*x[j];
        x[i] = (b[i]-sum)/a[i][i];
    }
    printf("线性方程组的解为: \n");
    for(int i = 1; i <= n; i++){
        if(x[i] == 0)x[i] = fabs(x[i]);
        printf("x[%d]=%f\n",i,x[i]);
    }
    return 0;
}

测试数据
10
4 2 -3 -1 2 1 0 0 0 0 5
8 6 -5 -3 6 5 0 1 0 0 12
4 2 -2 -1 3 2 -1 0 3 1 3
0 -2 1 5 -1 3 -1 1 9 4 2
-4 2 6 -1 6 7 -3 3 2 3 3
8 6 -8 5 7 17 2 6 -3 5 46
0 2 -1 3 -4 2 5 3 0 1 13
16 10 -11 -9 17 34 2 -1 2 2 38
4 6 2 -7 13 9 2 0 12 4 19
0 0 -1 8 -3 -24 -8 6 3 -1 -21
输出结果
这里写图片描述

列主元素高斯消元

列主元素消元是对高斯消元的该进,小主元可能导致计算失败,因为这个绝对值很小的数用作除数,乘数很大,引起约化中间结果数量级严重增长,所以在待消元的所在列中每一次选取最大的作为消元的主元素,并且交换第k行和第maxi行,再继续消元

#include <stdio.h> 
#include <stdlib.h>
#include <cmath> 
#define N 1000
double a[N][N],b[N],x[N];

int main(){
    FILE *in;
    int n;
    if((in = fopen("in.txt","r")) == NULL){printf("cannot open the file\n");exit(0);}
    fscanf(in,"%d",&n);
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n; j++)
            fscanf(in,"%lf",&a[i][j]);
        fscanf(in,"%lf",&b[i]);
    }
    bool flag = true;
    for(int k = 1; k < n; k++){
        int maxi;
        double ans = a[k][k],temp;
        for(int i = k+1; i <= n; i++)if(fabs(a[i][k]) > ans){ans = fabs(a[i][k]);maxi = i;} //在待消元的所在列中每一次选取最大的作为消元的主元素
        for(int j = k; j <= n; j++){temp = a[k][j]; a[k][j] = a[maxi][j]; a[maxi][j] = temp;}  //交换第k行和第maxi行
        temp = b[k],b[k] = b[maxi],b[maxi] = temp;
        if(a[k][k] == 0){printf("gauss method does not run"); flag = false; break;}  //如果a[k][k]为0 ,则不满足消元条件
        else {
            for(int i = k+1; i <= n; i++){
                for(int j = k+1; j <= n; j++)
                    a[i][j] = a[i][j]-a[k][j]*a[i][k]/a[k][k];
                b[i] = b[i]-b[k]*a[i][k]/a[k][k];
            }
        }
    }
    if(flag){
        for(int i = n; i > 0; i--){   //注意:  回代的时候要求a[n][n] != 0
            double sum=0;
            for(int j = i+1; j <= n; j++) sum = sum+a[i][j]*x[j];
            x[i] = (b[i]-sum)/a[i][i];
        }
    }
    printf("线性方程组的解为: \n");
    for(int i = 1; i <= n; i++)printf("x[%d] = %f\n",i,x[i]);
    return 0;
}

普通平方根法求解对称正定矩阵

对于一些特殊的矩阵,例如对称正定矩阵,平方根法比较实用(又叫Cholesky法)

#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <string.h>
#include <cmath>

using namespace std;
const int maxn = 100+10;
double a[maxn][maxn],b[maxn],l[maxn][maxn],x[maxn],y[maxn];

int main(){
    FILE *in;
    int n;
    if((in = fopen("in.txt","r")) == NULL){cout<<"cannot open the file"<<endl;exit(0);}
    fscanf(in,"%d",&n);
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n; j++){
            fscanf(in,"%lf",&a[i][j]);
        }
    }
    for(int i = 1; i <= n; i++)fscanf(in,"%lf",&b[i]);
    double sum;
    for(int i = 1; i <= n; i++){   //平方根法求L*L(T)
        for(int j = 1; j <= i-1; j++){  // i = j+1, j+2 ...
            sum = 0;
            for(int k = 1; k <= j - 1; k++) sum += l[i][k] * l[j][k];//一定要注意这里不是l[i][k]*l[k][j]
            l[i][j] = (a[i][j] - sum)/l[j][j];
        }
        sum = 0;
        for(int k = 1; k <= i-1; k++)sum += (l[i][k]*l[i][k]);
        l[i][i] = sqrt(a[i][i] - sum);
    }
    cout<<"L矩阵如下 :"<<endl;
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n; j++){
            cout<<l[i][j]<<" ";
        }
        cout<<endl;
    }
    for(int i = 1; i <= n; i++){
        double sum = 0;
        for(int j = 1; j <= i-1; j++)sum += l[i][j]*y[j];
        y[i] = (b[i] - sum)/l[i][i];
    }
    for(int i = 1; i <= n; i++)cout<<y[i]<<" ";cout<<endl;
    for(int i = 1; i <= n; i++){   //求转置矩阵
        for(int j = i+1; j <= n; j++){
            double temp = l[i][j];
            l[i][j] = l[j][i];
            l[j][i] = temp;
        }
    }
    for(int i = n; i >= 1; i--){
       double sum = 0;
        for(int j = i+1; j <= n; j++)sum += l[i][j]*x[j];
        x[i] = (y[i] - sum)/l[i][i];
    }
    for(int i = 1; i <= n; i++)cout<<x[i]<<" ";
    cout<<endl;
    return 0;
}

改进的平方根求法解对称正定矩阵

思想就是分解成A = LDL(T) ,其中D是对角阵,L是单位下三角阵,然后求解

#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <string.h>
using namespace std;

const int maxn = 4;
int main(){
    FILE *in;
    int n;
    double a[maxn][maxn],b[maxn],d[maxn],l[maxn][maxn],x[maxn],y[maxn];
    double sum;
    if((in = fopen("in.txt","r")) == NULL){cout<<"cannot open the file"<<endl;exit(0);}
    fscanf(in,"%d",&n);
    memset(l,0,sizeof(l));
    memset(d,0,sizeof(d));
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n; j++)
            fscanf(in,"%lf",&a[i][j]);
    }
    for(int i = 1; i <= n; i++)fscanf(in,"%lf",&b[i]);
    fclose(in);
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= i-1; j++){
            double sum = 0;
            for(int k = 1; k <= j - 1; k++) sum += d[k]*l[i][k]*l[j][k];
            l[i][j] = (a[i][j] - sum)/d[j];
        }
        sum = 0;
        for(int k = 1; k <= i-1; k++) sum += d[k]*l[i][k]*l[i][k];
        d[i] = a[i][i]-sum;
    }
    cout<<"D矩阵是:"<<endl;
    for(int i = 1; i <= n; i++)cout<<d[i]<<" ";//只有主对角线上有元素的对角阵
    cout<<endl;
    for(int i = 1; i <= n; i++)l[i][i] = 1;  //把对角线上的赋值为1
    cout<<"L矩阵是:"<<endl;
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n; j++)
            cout<<l[i][j]<<" ";
        cout<<endl;
    }
    for(int i = 1; i <= n; i++){
        sum = 0;
        for(int k = 1; k <= i-1; k++) sum += l[i][k]*y[k];
        y[i] = b[i]-sum;
    }
    for(int i = n; i >= 1; i--){
        sum = 0;
        for(int k = i+1; k <= n; k++) sum += l[k][i]*x[k];
        x[i] = y[i]/d[i]-sum;
    }
    for(int i = 1; i <= n; i++)
        printf("x[%d] = %lf\n",i,x[i]);
    return 0;
}

对称正定矩阵测试例子
3
1 1 2
1 2 0
2 0 11
5 8 7
输出
这里写图片描述

矩阵三角分解(LU)(杜立特分解)

这个也是比较经典的算法,定理:如果A的各阶主子式不为0,A可以唯一的分解成一个单位下三角L和一个非奇异的上三角U的乘积,分解之后求解

#include <iostream>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>

using namespace std;
const int maxn = 10+10;

double a[maxn][maxn],L[maxn][maxn],U[maxn][maxn];
double b[maxn],x[maxn],y[maxn];

int main(){
    int n;
    //freopen("in.txt","r",stdin);
    cin>>n;
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n; j++){
            cin>>a[i][j];
        }
    }
    for(int i = 1; i <= n; i++)cin>>b[i];
    memset(L,0,sizeof(L));
    memset(U,0,sizeof(U));
    for(int i = 1; i <= n; i++)U[1][i] = a[1][i];
    for(int i = 1; i <= n; i++)L[i][1] = a[i][1]/U[1][1];
    for(int k = 1; k <= n; k++){
        for(int j = k; j <= n; j++){    //用L的第K行去乘U的每一列,得到U的第K行,j是第K行的第j列
            double sum = 0;
            for(int i = 1; i < k; i++)   //i是用来计算第j列之前的数的和
                sum = sum+L[k][i]*U[i][j];
            U[k][j] = a[k][j] - sum;
        }
        for(int i = k; i <= n; i++){ //用L的每一行去乘U的第K列,得到L的第K列,其中i是第K列的第i行
            double sum = 0;
            for(int j = 1; j < k; j++) //j是指第i列前面的数的和
                sum = sum + L[i][j]*U[j][k];
            L[i][k] = (a[i][k]-sum)/U[k][k];
        }
    }
    cout<<"L矩阵如下:"<<endl;
    for(int i = 1; i <= n; i++){//L矩阵
        for(int j = 1; j <= n; j++){
            cout<<L[i][j]<<" ";
        }
        cout<<endl;
    }
    cout<<"U矩阵如下:"<<endl;
    for(int i = 1; i <= n; i++){ //U矩阵
        for(int j = 1; j <= n; j++){
            cout<<U[i][j]<<" ";
        }
        cout<<endl;
    }
     y[1] = b[1];
     for(int i = 2; i <= n; i++){
        double sum = 0;
        for(int k = 1; k < i; k++)sum += L[i][k]*y[k];
        y[i] = b[i] - sum;
     }
     x[n] = y[n]/U[n][n];
     for(int i = n - 1; i >= 1; i--){
        double sum = 0;
        for(int k = i+1; k <= n; k++)sum += U[i][k]*x[k];
        x[i] = (y[i]-sum)/U[i][i];
     }
     cout<<"线性方程组的解是:"<<endl;
     for(int i = 1; i <= n; i++)cout<<x[i]<<" ";
     cout<<endl;
}

测试数据
3
1 2 3
1 3 5
1 3 6
2 4 5
输出
这里写图片描述


追赶法 求解 三对角方程组

求解三对角方程组方程

#include <iostream>
#include <stdio.h>
const int maxn = 10 + 10;
using namespace std;

double a[maxn], b[maxn], c[maxn],f[maxn];
double l[maxn], u[maxn];
double x[maxn], y[maxn];

int main() {
    freopen("in.txt", "r", stdin);
    freopen("out.txt", "w", stdout);
    int n;
    scanf("%d", &n); 
    for (int i = 1; i <= n - 1; i++)scanf("%lf", &c[i]);//三对角的上方数组
    for (int i = 2; i <= n; i++)scanf("%lf", &a[i]);//三对角的下方数组
    for (int i = 1; i <= n; i++)scanf("%lf", &b[i]);//三对角的中间数组
    for (int i = 1; i <= n; i++)scanf("%lf", &f[i]);//右边的矩阵
    l[1] = b[1];
    for (int i = 1; i <= n - 1; i++) {//求解得到l数组和u数组
        u[i] = c[i] / l[i];
        l[i + 1] = b[i + 1] - a[i + 1] * u[i];
    }
    y[1] = f[1] / l[1];
    for (int i = 2; i <= n; i++)y[i] = (f[i] - a[i] * y[i - 1]) / l[i];//分两步,先求解y数组
    x[n] = y[n];
    for (int i = n - 1; i >= 1; i--)x[i] = y[i] - u[i] * x[i + 1];//然后求解x数组
    for (int i = 1; i <= n; i++)printf("x[%d] = %6.4f\n", i, x[i]);
    return 0;
}

输入输出
/*
输入:
4
-1 -2 -2
-1 -2 -2
2 3 4 5
3 1 0 -5
*/
这里写图片描述
/*
输出:
x[1] = 2.0000
x[2] = 1.0000
x[3] = 0.0000
x[4] = -1.0000
*/


雅可比迭代法

雅可比求解线程方程组和之前的不同,是一种迭代法来求解线性方程组的解

#include <iostream>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>

using namespace std;
const int M = 100;  //最大的迭代次数
const int maxn = 100;

double a[maxn][maxn],b[maxn],x[maxn],y[maxn];
int n;

int main(){
    freopen("in.txt","r",stdin);
    memset(x,0,sizeof(x));//选取初值全为0
    memset(y,0,sizeof(y));
    cin>>n;
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n; j++){
            cin>>a[i][j];
        }
    }
    for(int i = 1; i <= n; i++)cin>>b[i];
    double sum, maxx, eps = 0.000000001;
    bool flag = true;
    int k;
    for(k = 1; k <= M; k++){  //控制迭代次数
        for (int i = 1; i <= n ; i++ ){ //对i = 1,2,3,..n计算公式
            sum = 0;
            for (int j = 1; j <= n; j++)
                if(j != i) sum += a[i][j] * x[j];
            y[i] = ( b[i] - sum) / a[i][i]; //由x(k)计算出x(k+1)
        }
        maxx = 0;
        for(int i = 1; i <= n; i++)
            if( maxx < fabs(x[i]-y[i])) maxx = fabs(x[i] - y[i]);//当相邻两次迭代结果的偏差max(x[i] - y[i]) < eps即小于给定的精度可以输出
        if ( maxx < eps ){flag = false; break;}
        printf( "\nk = %d, ",k);
        for(int i = 1; i <= n; i++ ) printf( "y[%d] = %lf ",i,y[i]);
        for(int i = 1; i <= n; i++ ) x[i] = y[i];  //否则继续迭代
    }
    if(flag)printf( "ERROR!\n");
    else {
        printf( "\nk = %d, ", k);
        for (int i = 1; i <= n; i++ ) printf( "y[%d] = %lf ", i, y[i] );
    }
    return 0;
}

测试数据
3
8 -3 2
4 11 -1
6 3 12
20 33 36
输出
这里写图片描述

高斯赛德尔迭代法及超松弛迭代法

这两个方法是雅可比迭代法的改进

#include <iostream>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>

using namespace std;
const int M = 100;  //最大的迭代次数
const int maxn = 100;

double a[maxn][maxn],b[maxn],x[maxn],y[maxn];
int n;

int main(){
    freopen("in.txt","r",stdin);
    memset(x,0,sizeof(x));//选取初值全为0
    memset(y,0,sizeof(y));
    cin>>n;
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n; j++){
            cin>>a[i][j];
        }
    }
    for(int i = 1; i <= n; i++)cin>>b[i];
    double sum, maxx, eps = 0.000000001,w = 1;//当w = 1时就是高斯赛德尔迭代法
    bool flag = true;
    int k;
    for(k = 1; k <= M; k++){  //控制迭代次数
        for (int i = 1; i <= n ; i++ ){ //对i = 1,2,3,..n计算公式
            sum = 0;
            for(int i = 1; i <= n; i++ ) y[i] = x[i];//先把原来的值存好,等下判断退出的时候要用
            for (int j = 1; j <= n; j++)
                if(j != i) sum += a[i][j] * x[j];
            x[i] = ( b[i] - sum) / a[i][i]; //由x(k)计算出x(k+1)
            x[i] = (1 - w)*y[i]+w*x[i];//y[i]存的是原来的值
        }
        maxx = 0;
        for(int i = 1; i <= n; i++)
            if( maxx < fabs(x[i]-y[i])) maxx = fabs(x[i] - y[i]);//当相邻两次迭代结果的偏差max(x[i] - y[i]) < eps即小于给定的精度可以输出
        if ( maxx < eps ){flag = false; break;}
        printf( "\nk = %d, ",k);
        for(int i = 1; i <= n; i++ ) printf( "y[%d] = %lf ",i,y[i]);
    }
    if(flag)printf( "ERROR!\n");
    else {
        printf( "\nk = %d, ", k);
        for (int i = 1; i <= n; i++ ) printf( "y[%d] = %lf ", i, y[i] );
    }
    return 0;
}

拉格朗日插值

这个是根据一些已知点求解近似函数的算法
这里写图片描述

#include <iostream>
#include <stdio.h>
#include <string.h>
using namespace std;

const int maxn = 100 + 10;
int main(){
    //freopen("in.txt","r",stdin);
    int n;
    double x[maxn],y[maxn],X;
    scanf("%d",&n);
    for(int i = 0; i < n; i++)scanf("%lf",&x[i]);
    for(int i = 0; i < n; i++)scanf("%lf",&y[i]);
    while(scanf("%lf",&X) != EOF){
        double sum = 0;
        for(int i = 0; i < n; i++){
            double temp = 1;
            for(int j = 0; j < n; j++){
                if(i != j){
                    temp *= (X - x[j])/(x[i] -x[j]);
                }
            }
            sum += temp*y[i];
        }
        printf("%lf\n",sum);
    }
    return 0;
}
/*
6
0.4 0.55 0.65 0.80 0.95 1.05
0.41075 0.57815 0.69675 0.90 1.00 1.25382
0.596
*/

最小二乘法求解拟合曲线

这个也是求解拟合曲线的另一种方法,过程看一下书吧,直接按书上的思路撸的代码

//最小二乘法求解拟合曲线
//求解过程看书上P116(刘师少版)

#include <iostream>
#include <stdio.h>
#include <math.h>
using namespace std;
const int maxn = 100 + 10;

double x[maxn], y[maxn];//数据的
double a[maxn][maxn], b[maxn];//高斯消元中的矩阵
double ans[maxn];//要求解的a数组
double R[maxn];

int main() {
    //freopen("in.txt", "r", stdin);
    //freopen("out.txt", "w", stdout);
    int n, m;
    scanf("%d", &n);//数据的组数
    for (int i = 1; i <= n; i++)scanf("%lf", &x[i]);
    for (int i = 1; i <= n; i++)scanf("%lf", &y[i]);
    scanf("%d", &m);//多项式拟合的次数  : 例如 1 代表 直线拟合  a0 + a1*x 
    a[0][0] = n;
    for (int i = 1; i <= m; i++) {//计算增广矩阵第一行
        double sum = 0;
        for (int j = 1; j <= n; j++)sum += pow(x[j], i);
        a[0][i] = sum;
    }
    for (int i = 1; i <= m; i++) {
        for (int j = 0; j <= m; j++) {
            double sum = 0;
            for (int k = 1; k <= n; k++)sum += pow(x[k], j + i);
            a[i][j] = sum;
        }
    }
    for (int i = 0; i <= m; i++) {
        double sum = 0;
        for (int j = 1; j <= n; j++)sum += y[j] * pow(x[j], i);
        b[i] = sum;
    }
    int k;//高斯消元求解
    for (k = 0; k < m; k++) {//m - 1次消元
        for (int i = k + 1; i <= m; i++) {
            for (int j = k + 1; j <= m; j++) ///注意一定要从 K + 1开始
                a[i][j] = a[i][j] - a[k][j] * a[i][k] / a[k][k];
            b[i] = b[i] - b[k] * a[i][k] / a[k][k];
        }
    }
    for (int i = m; i >= 0; i--) {//求解ans数组(逆向)
        double sum = 0;
        for (int j = i + 1; j <= m; j++) sum = sum + a[i][j] * ans[j];
        ans[i] = (b[i] - sum) / a[i][i];
    }
    printf("a数组如下:\n");
    for (int i = 0; i <= m; i++) printf("a[%d] = %6.4f\n", i, ans[i]);
    for (int i = 1; i <= n; i++) {
        double sum = ans[0];
        for (int j = 1; j <= m; j++)sum += ans[j] * pow(x[i], j);
        R[i] = sum - y[i];
    }
    printf("\n余项R数组为:\n");
    for (int i = 1; i <= n; i++)printf("R[%d] = %6.4f\n", i, R[i]);
    return 0;
}
/*
12
0 5 10 15 20 25 30 35 40 45 50 55
0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64
3
*/
/*
4
1.36 1.73 1.95 2.28
14.094 16.844 18.475 20.963
1
输出:
a数组如下:
a[0] = 3.9374
a[1] = 7.4626
余项R数组为:
R[1] = -0.0074
R[2] = 0.0037
R[3] = 0.0145
R[4] = -0.0108
*/
/*
6
0 1 2 3 4 5
5 2 1 1 2 3
2
输出:
a数组如下:
a[0] = 4.7143
a[1] = -2.7857
a[2] = 0.5000
余项R数组为:
R[1] = -0.2857
R[2] = 0.4286
R[3] = 0.1429
R[4] = -0.1429
R[5] = -0.4286
R[6] = 0.2857
*/

牛顿插值多项式 (均差)

这个是比较经典的,先求均差,然后求解差值函数,代码如下:

//牛顿插值多项式求解插值函数
//公式以及解释见书上 P100,P101
#include <iostream>
#include <stdio.h>
#include <string.h>
#include <string>

using namespace std;
const int maxn = 100 + 10;
double x[maxn], y[maxn], F[maxn][maxn],X;

int main() {
    //freopen("in.txt", "r", stdin);
    //freopen("out.txt", "w", stdout);
    int n, m;
    scanf("%d", &n);
    for (int i = 0; i < n; i++)scanf("%lf", &x[i]);
    for (int i = 0; i < n; i++)scanf("%lf", &y[i]);
    for (int i = 0; i < n; i++)F[i][0] = y[i];//二维数组的第一列
    for (int j = 1; j < n; j++) {//求均差的公式
        for (int i = j; i < n; i++) {
            F[i][j] = (F[i][j - 1] - F[i - 1][j - 1]) / (x[i] - x[i - j]);
        }
    }
    printf("均差表如下:\n");
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            printf("%6.4f", F[i][j]);
        }
        printf("\n");
    }
    while (scanf("%lf", &X) != EOF) {//输入要求的x对应的函数值
        double sum = 0;
        for (int i = 0; i < n; i++) {//求函数值
            double temp = F[i][i];
            for (int j = 0; j < i; j++)temp *= (X - x[j]);
            sum += temp;
        }
        printf("\n%8.6f\n", sum);
    }
    return 0;
}
/*
输入:
5
0.4 0.55 0.65 0.8 0.9
0.4175 0.57815 0.69657 0.88811 1.02652
0.5
0.85
1.05
*/
/*
输出:
均差表如下:
0.41750.00000.00000.00000.0000
0.57821.07100.00000.00000.0000
0.69661.18420.45280.00000.0000
0.88811.27690.3709-0.20470.0000
1.02651.38410.42870.16500.7392
0.522016
0.956050
1.258229
*/

变步长梯形求积法

这个是一个求积分的算法,基本思想:在求积的工程中,通过对计算结果精度的不断估计,逐步改变步长(逐次分半)直到精度满足要求为止,求解过程也不细说,看书。

//刘师少书本P138
#include <iostream>
#include <stdio.h>
#include <math.h>
using namespace std;
const int maxn = 100 + 10;
double a,b,n,eps = 0.00001;

double f(double x){
    //return (x/(4+x*x));
    return 4/(1+x*x);
    //return sin(x);
    //return sin(x)/x;
}

int main(){
    freopen("in.txt","r",stdin);
    while(scanf("%lf%lf",&a,&b) != EOF){  //积分上下限以及等分数
        double h = b - a;   //一开始的等分区间(分成一块)
        double T1 = (h/2)*(f(a) + f(b)),T2;   //最初的梯形求积公式
        int n = 1;
        for(;;){
            double S = 0;
            for(double x = a + h/2; x < b; x += h)S += f(x);
            T2 = (T1 + h*S)/2;
            printf("步长为 %d 时, Tn = %8.6f, T2n = %8.6f\n",n,T1,T2);
            if(fabs(T2 - T1) < eps)break;
            T1 = T2,h /= 2,n *= 2;//n记录步长
        }
        printf("\n经过变步长梯形求积公式得到的积分结果是 :%8.6f\n",T2);
    }
    return 0;
}
/*
输入
0 1
*/
/*
输出
0.111566
*/

复化梯形公式 与 复化辛普森公式求积

这两个也是比较好的求积公式,具体看书

//复化梯形求积法的算法实现
//计算方法刘师少P134(梯形公式),P135(辛普森)
#include <iostream>
#include <stdio.h>
#include <math.h>
using namespace std;
const int maxn = 100 + 10;
double a,b,n;

double f(double x){
    return sin(x);
}

double FuHuaTiXing(double a,double b,double n){//复化梯形求积公式
    double h = (b - a)/n;//步长
    double T = 0;
    for(int i = 1; i <= n - 1; i++)T += f(a + i*h);
    T = (h/2)*(f(a) + 2*T + f(b));
    return T;
}

double FuHuaXingPuSeng(double a,double b,double n){//复化辛普森求积公式
    double h = (b - a)/n;
    double S1 = f(a + h/2),S2 = 0,S;
    for(int i = 1; i <= n - 1; i++){
        S1 += f(a + i*h + h/2);
        S2 += f(a + i*h);
    }
    S = (h/6)*(f(a) + 4*S1 + 2*S2 + f(b));
    return S;
}

int main(){
    freopen("in.txt","r",stdin);
    while(scanf("%lf%lf%lf",&a,&b,&n) != EOF){//积分上下限以及等分数
        printf("复化梯形求出来的积分是: %6.6f\n",FuHuaTiXing(a,b,n));
        printf("复化辛普森求出来的积分是: %6.6f\n",FuHuaXingPuSeng(a,b,n));
    }
    return 0;
}
/*
输入
1 2 8
*/
/*
输出
0.955203
0.956449
*/

改进的欧拉法 求解 微分方程的 数值

欧拉法求解微分方程

//  常微分方程的数值解法
//  改进的欧拉公式,具体的解释看 刘师少版计算方法 P153 - P156
//  y[i + 1]' = yi + h * f(xi,yi);//预测(欧拉公式)
//  y[i + 1] = yi + h * (f(xi,yi) + f(x[i + 1],y[i+1]'))// (梯形公式校正)

#include <iostream>
#include <stdio.h>
#include <string.h>

using namespace std;
const int maxn = 100 + 10;

double f(double x, double y) {//P158例题
    return y - 2 * x / y;
}

int main() {
    freopen("in.txt", "r", stdin);
    //freopen("out.txt", "w", stdout);
    double x0, y0, xn, x1, y1;
    //double yp, yc;
    int n;
    scanf("%lf%lf%lf%d", &x0, &y0, &xn, &n);//输入区间左边得x0,y0,和区间右边的xn,以及区间数
    double h = (xn - x0) / n;  //求步长
    for (int i = 1; i <= n; i++) {
        x1 = x0 + h; //求出下一个点的 x 值
        //yp = y0 + h * f(x0, y0);
        //yc = y0 + h * f(x1, yp);
        //y1 = (yp + yc) / 2;
        y1 = y0 + h / 2 * (f(x0, y0) + f(x1,y0 + h * f(x0, y0)));
        printf("x[%d] = 1%8.6f   y[x[%d]] = %8.6f\n",i,x1,i, y1);
        x0 = x1;
        y0 = y1;
    }
    return 0;
}

/*
输入:
0 1 1 10
*/
/*
输出:
x[1] = 10.100000   y[x[1]] = 1.095909
x[2] = 10.200000   y[x[2]] = 1.184097
x[3] = 10.300000   y[x[3]] = 1.266201
x[4] = 10.400000   y[x[4]] = 1.343360
x[5] = 10.500000   y[x[5]] = 1.416402
x[6] = 10.600000   y[x[6]] = 1.485956
x[7] = 10.700000   y[x[7]] = 1.552514
x[8] = 10.800000   y[x[8]] = 1.616475
x[9] = 10.900000   y[x[9]] = 1.678166
x[10] = 11.000000   y[x[10]] = 1.737867
*/

数值积分 之 龙贝格求积分法

这个也是求解数值积分的一个好算法,具体过程看书

//龙贝格求积法求解积分
//f函数来自书上的例题,计算结果正确
#include <iostream>
#include <stdio.h>
#include <math.h>
#include <string.h>

using namespace std;
const int maxn = 100 + 10;
double a,b,eps = 1e-5;
double num[maxn][4];
int k;

double f(double x){
    return 4/(1+x*x);
}

void Romberg(double a,double b){
    k = 1;
    double T1,T2,S1,S2,C1,C2,R1,R2;
    double h = b - a;//一开始的区间
    T1 = (h/2)*(f(b) + f(a));
    num[0][0] = T1;
    for(;;){
        double sum = 0;
        for(double x = a + h/2; x < b; x += h)sum += f(x);
        T2 = (T1 + sum*h)/2;
        num[k][0] = T2;
        S2 = T2 + (T2 - T1)/3;//梯形
        num[k][1] = S2;
        if(k != 1){
             C2 = S2 + (S2 - S1)/15;//辛普森
             num[k][2] = C2;
             if(k != 2){
                R2 = C2 + (C2 - C1)/63; //龙贝格
                num[k][3] = R2;
                if(k != 3)if(fabs(R2 - R1) < eps){printf("龙贝格求积求得结果是:  %8.6f\n\n",R2);break;}
                R1 = R2;
            }
            C1 = C2;
        }
        h /= 2;k++;//逐步减小步长
        T1 = T2;
        S1 = S2;
    }
}

int main(){
    //freopen("in.txt","r",stdin);
    //freopen("out.txt","w",stdout);
    while(scanf("%lf%lf",&a,&b) != EOF){//输入积分上下限
        memset(num,0,sizeof(num));
        Romberg(a,b);
        printf("龙贝格算法计算表格如下:\n");
        printf("%c%6c%12c%12c%12c\n",'k','T','S','C','R');
        for(int i = 0; i < k; i++){
            printf("%d  ",i);
            for(int j = 0; j < 4; j++){
                if(num[i][j] != 0)
                printf("%8.6f    ",num[i][j]);
            }
            printf("\n");
        }
    }
    return 0;
}
/*
输入
0 1
*/
/*
输出
龙贝格求积求得结果是:  3.141593
龙贝格算法计算表格如下:
k     T           S           C           R
0  3.000000
1  3.100000    3.133333
2  3.131176    3.141569    3.142118
3  3.138988    3.141593    3.141594    3.141586
*/

常微分方程初值问题 数值解法 之 经典 - 四阶龙格库塔 法

常微分方程数值解法,这个也是比较经典的

//经典的四阶龙格 - 库塔 法 求解常微分方程的数值解
//算法的具体解释 等件  刘师少版计算方法 P159 - P163
#include <iostream>
#include <stdio.h>
#include <string.h>

using namespace std;
const int maxn = 100 + 10;

double f(double x, double y) {// P162 的 例题
    return 2 * x * y;
}

int main() {
    //freopen("in.txt", "r", stdin);
    //freopen("out.txt", "w", stdout);
    double x0, y0, xn,x1,y1;
    double k1, k2, k3, k4;//四个中间的斜率
    int n;
    scanf("%lf%lf%lf%d", &x0, &y0, &xn, &n);//输入左边区间的  x0 ,y0,以及右边区间的 xn值,以及分的区间的个数
    double h = (xn - x0) / n;  //步长
    for (int i = 1; i <= n; i++) {
        x1 = x0 + h;
        k1 = f(x0, y0);
        k2 = f(x0 + h / 2, y0 + h * k1 / 2);
        k3 = f(x0 + h / 2, y0 + h * k2 / 2);
        k4 = f(x1, y0 + h * k3);
        y1 = h * (k1 + 2 * k2 + 2 * k3 + k4)/6 + y0;
        printf("x[%d] = %8.6f  y[%d] = %8.6f\n", i, x1, i, y1);
        x0 = x1, y0 = y1;
    }
    return 0;
}
/*
输入:
0 1 1 5
*/
/*
输出:
x[1] = 0.200000  y[1] = 1.040811
x[2] = 0.400000  y[2] = 1.173510
x[3] = 0.600000  y[3] = 1.433321
x[4] = 0.800000  y[4] = 1.896441
x[5] = 1.000000  y[5] = 2.718107
*/

阅读更多
想对作者说点什么? 我来说一句

没有更多推荐了,返回首页