目录

• 矩阵相乘
• 各种求解一元线性方程解的解法
• 普通高斯消元
• 列主元素高斯消元
• 普通平方根法求解对称正定矩阵
• 改进的平方根求法解对称正定矩阵
• 矩阵三角分解（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; //倒数为0double 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] != 0double 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;  //把对角线上的赋值为1cout<<"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));//选取初值全为0memset(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));//选取初值全为0memset(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
*/

---

常用计算方法(C语言代码)(计算方法课程)相关推荐

1. C语言编程计算差商表,计算方法C语言编程计算方法C语言编程.doc

   计算方法C语言编程第二章2已知方程在区间[1,2]内有一根,试问用二分法求根,使其具有5位有效数字至少应二分多少次?[程序设计]#includemain(){int n=0; float x1=1.0 ...

2. 计算方法c语言编程,计算方法C语言编程计算方法C语言编程.doc

   计算方法C语言编程第二章2已知方程在区间[1,2]内有一根,试问用二分法求根,使其具有5位有效数字至少应二分多少次?[程序设计]#includemain(){int n=0; float x1=1.0 ...

3. c语言math函数 sgn,常用矩阵计算C语言代码

   参考资料: 行列式:http://zh.wikipedia.org/wiki/行列式#.E4.BB.A3.E6.95.B0.E4.BD.99.E5.AD.90.E5.BC.8F 伴随矩阵:http:/ ...

4. 用c语言编程计算10,计算方法c语言编程.doc

   计算方法c语言编程 计算方法C语言编程 1.已知方程在区间[1,2]内有一根,试问用二分法求根,使其具有5位有效数字至少应二分多少次? [程序设计] clc;clear a=1;b=2;n=0; x= ...

5. 计算方法(c语言版)靳天飞答案,《数值分析教程》PPT课件.ppt

   <<数值分析教程>PPT课件.ppt>由会员分享,可在线阅读,更多相关<<数值分析教程>PPT课件.ppt(47页珍藏版)>请在装配图网上搜索. 1.计 ...

6. 图形基本变换c语言代码,图形变换-C语言课程设计.doc

   学号 <> 课程设计报告 图形变换网络工程班级:16(3)姓名:指导教师:成绩: 计算机学院 2017 年 5月 10日 目录- 1 - 1 设计要求- 2 - 2 程序功能- 2 - 3 ...

7. 等额本息还款方式的年利率计算方法及java代码实现

   等额本息还款方式在生活中是很常用的,比如贷款买房,买车,信用卡分期还款,京东白条分期 ,等等,只要是借了钱,每月还款金额一样的,都属于等额本息的方式.有些时候我们贷款后,商家或银行告诉了每期的还款额, ...

8. 浅析VS2010反汇编 VS 反汇编方法及常用汇编指令介绍 VS2015使用技巧 调试-反汇编 查看C语言代码对应的汇编代码...

   浅析VS2010反汇编 2015年07月25日 21:53:11 阅读数:4374 第一篇 1. 如何进行反汇编 在调试的环境下,我们可以很方便地通过反汇编窗口查看程序生成的反汇编信息.如下图所示. ...

9. c语言编程计算c上0下n,计算方法C语言编程讲解.doc

   计算方法C语言编程讲解 计算方法C语言编程 1.已知方程在区间[1,2]内有一根,试问用二分法求根,使其具有5位有效数字至少应二分多少次? [程序设计] clc;clear a=1;b=2;n=0; ...

最新文章

1. js 调用C#.NET后台方法 转载自：http://www.cnblogs.com/lizhao/archive/2010/11/23/1990436.html...
2. [原创] Android SDK 安装全记录
3. c++中静态成员变量和静态成员函数
4. oracle数据库tx锁,oracle数据库有把TX锁，如何定位锁在哪？
5. 软件项目管理与素质拓展-2.2什么是项目
6. C# WinForm窗体程序、如何实现像QQ一样的热键
7. SQLiteSpy介绍和使用
8. win7计算机无法连接投影仪,笔记本win7系统连接投影仪显示没信号如何解决？
9. python工程师认证证书报考条件_Python工程师需要具备什么条件
10. WireShark 抓包使用教程--详细
11. CPU处理器的分类（ARM系列中央处理器）
12. 如何看待阿里云推出的免费虚拟主机？
13. 第二届国际计算成像会议-CITA
14. OC textField键盘弹起事件
15. vue前端项目启动出错处理
16. 360手机安全专家给出了八大安全防范建议
17. 开发人员各级岗位胜任力模型
18. 三大通识知识：进程，线程，网络(四)
19. android10 保存图片到系统相册，刷新媒体库
20. 宽带连接远程计算机没反应6,宽带连接正在连接通过WAN微型端口错误678的解决办法...

热门文章

1. PClint 使用教程
2. kettle入门(二) 之 kettle连接oracle报的坑爹错误 Error occured while trying to connect to the database 的几种情况
3. 黑客劫持域名步骤大曝光
4. neovim初始化以及插件安装
5. jsp实现数据提交以及jsp数据保存到本地
6. 大数据架构师之路-大数据框架大全
7. web页面播放优酷视频，播放html5视频，兼容ie7 vcastr22.swf播放
8. golang学习笔记之string转换
9. Android studio修改标题菜单栏增加功能图标（navigation bar toolbar）
10. ratingbar 的使用