目录
- 矩阵相乘
- 各种求解一元线性方程解的解法
- 普通高斯消元
- 列主元素高斯消元
- 普通平方根法求解对称正定矩阵
- 改进的平方根求法解对称正定矩阵
- 矩阵三角分解(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
*/