数值分析

拉格朗日插值#include<iostream>using namespace std;unsigned int i,j,k,n;float a[200]={0},b[200]={0};//数组a用来存放x[i]b用来存放y[i]double x,xmt=1,sum=0;//xmt为分式累乘器，sum为累加器void input()//输入函数{cout<<"Please input n"<<endl<<"n=";cin>>n;if(n>199)//防止数组溢出{cout<<"unable to calculate,n is too big"<<endl<<"Please type in n"<<endl<<"n=";cin>>n;}else{for(i=0;i<=n;i++)//输入x[i]{cout<<"Please input x"<<i<<endl<<"x"<<i<<"=";cin>>a[i];}for(i=0;i<=n;i++)//输入y[i]{cout<<"Please input y"<<i<<endl<<"y"<<i<<"=";cin>>b[i];}cout<<"Please input x"<<endl<<"x=";//输入xcin>>x;}}double cal()//计算函数{for(k=0;k<=n;k++){for(j=0;j<=n;j++){if(j!=k){xmt*=(1.0*(x-a[j])/(a[k]-a[j]));}}sum+=b[k]*xmt;xmt=1;//累乘器清零}return(sum);}void main(){while(1)//让程序可以一直工作，而不是一次性的{input();cout<<"拉格朗日插值结果为:y="<<cal()<<endl;//结果sum=0;//累加器清零，便于下次计算}}牛顿插值#include<iostream>#include<math.h>void main（）{int n,i,j;double A[50][50];double x[50],y[50];double K=1,X=0,N=0,P;cout<<"请输入插值阶数："<;cin>>n;for (i=0;i<=n;i++){cout<<"请输入x"<<i<<"=";cin>>x[i];cout<<"请输入y"<<i<<"=";cin>>y[i];A[i][0]=x[i];A[i][1]=y[i];}for(j=2;j<=n+1;j++){for(i=1;i<=n;i++){A[i][j]=(A[i][j-1]-A[i-1][j-1])/(A[i][0]-A[i-j+1][0]);} }for(i=0;i<=n;i++){cout<<"输出第"<<i<<"阶插值为："<<A[i][i+1]<<endl;} cout<<"索要带入计算的x的值：X=";cin>>X;for(i=0;i<n;i++){K*=X-x[i];N+=A[i+1][i+2]*K;P+A[0][1]+N;}cout<<"所要求的函数值为：y="<<P<<endl;}龙贝格#include<iostream>#include<cmath>using namespace std;const int MAXRepeat = 100; //最大允许重复double function(double x)//被积函数,根据自己的需要手工输入{double s;s = 1.0 / (1 + x);return s;}void Romberg(double a, double b, double epsion, double f(double x)){int m = 1;int n = 1;int k;double h;double ep;double p;double xk;double s;double q;double T[MAXRepeat];h = b - a;T[0] = 0.5 * h * (f(a) + f(b));ep = epsion + 1.0;while ((ep >= epsion) && (m < MAXRepeat)){p = 0.0;for (k = 0; k < n; k++){xk = a + (k + 0.5) * h; // n-1p = p + f(xk); //计算∑f(xk+h/2),T} // k=0p = (T[0] + h * p) / 2.0; //T`m`(h/2),变步长梯形求积公式s = 1.0;for (k = 1; k <= m; k++){s = 4.0 * s; //[pow(4,m)T`m`(h/2)-T`m`(h)]/[pow(4,m)-1],2m阶牛顿柯斯特公式，即龙贝格公式q = (s * p - T[k - 1]) / (s - 1.0);T[k-1] = p;p = q;}ep = fabs(q - T[m - 1]);m++;T[m - 1] = q;n++; // 2 4 8 16h /= 2.0;}for (int i = 0; i < m; i++){int j;if (!(i % j)){cout<<T[i]<<endl;}else{cout<<T[i]<<" ";}j++;}}int main(){double a;double b;double epsion;cout<<"Please input the lower limit: ";cin>>a;cout<<"Please input the upper limit: ";cin>>b;cout<<"Please input the precision : ";cin>>epsion;Romberg( a, b, epsion, function);return 0;二分法#include <iostream.h>#include <math.h> #defineeps 0.1double fun(double x) { return x*x-3*x; }doubled iff_fun(double x) { return 2*x-3; }double dichotomy(double a,double b) {double c = 0.0;if((diff_fun(a) < 0) && (diff_fun(b) > 0)) { while(true) {c = (a + b)/2;if(diff_fun(c) < 0) {a=c;if(fabs((a-b)) <eps){ return (a+b)/2; }}else if(diff_fun(c) == 0){ return c; }else { b = c;if(fabs((a-b)) <eps){ return (a+b)/2; }}}}迭代法#include<iostream>#include<string>#include<cmath>using namespace std;double e=2.718281818284;double f(double x){double a=pow(e,-1*x);return x-(x-a)/(1+a);}void NewtonDiedai(double x,double d){ double a=x;double b=f(a);int k=0; //记录循环的次数while(((a-b)>d) || ((a-b)<-1*d)){ cout<<a<<endl;a=b;b=f(a);k++;if(k>100){cout<<"迭代失败！（可能是函数不收敛）"<<endl; return ;}}cout<<b<<endl;return;}int main(){cout<<"请输入初始值x0和要求得结果的精度：";double x,d;cin>>x>>d;NewtonDiedai(x,d);return 0;}/*雅可比算法的代码实现*/#include<iostream>#include<iomanip>#include<string>#include<vector>using namespace std;//函数求数组中的最大值double MaxOfList(vector<double>x){double max=x[0];int n=x.size();for(int i=0;i<n;i++)if(x[i]>max) max=x[i];return max;}迭代法的加速#include <iostream>#include<iomanip.h>using namespace std;int change_num=0;void display(double **l_array,int r_size1,int r_size2) //这个函数是用来显示矩阵（主要是用来调试使用的）{for(int i=0; i<r_size1; i++)for(int j=0; j<r_size2; j++){if(j!=0){cout<<setprecision(4)<<setiosflags(ios_base::left)<<setw(7)<<l_array[i][j];}else{cout<<setprecision(4)<<setiosflags(ios_base::left)<<l_array[i][j]<<" ";}if(j+1==r_size2){cout<<endl;}}cout<<resetiosflags(ios_base::left);}//这个下面的是算出上三角形bool Change(double **&array,int r_size,int c_size){int k=0;double max=array[k][k];int num=0;bool flag=false;while(k<r_size) //k表示现在在哪一行{max=array[k][k];for(int i=k; i<r_size; i++){if(array[i][k]>max){flag=true;max=array[i][k];num=i; //num表示最大行的行数}}if(!flag){num=k;}if(max==0) //如果一个列中的最大数是0的话那就说明这个矩阵的奇异矩阵{return false;}else if(k!=num){double temp=0;for(int i=0; i<c_size; i++){temp=array[k][i];array[k][i]=array[num][i];array[num][i]=temp;change_num++;}}for(int i=k; i<r_size-1; i++){double temp=array[i+1][k]; //这里的temp一定要保存的，因为直接使用array[i+1][k]，在下面的计算中会改变for(int j=k; j<c_size; j++){array[i+1][j]-=((temp/array[k][k])*array[k][j]);}}k++;flag=false;}return true;}double Rowlay(double **&array,int r_size, int c_size ) // 求行列式{double sum=1;for(int i=0; i<r_size; i++){sum*=array[i][i];}return sum;}void Equation(double **&array,int r_size, int c_size) //求方程组解{for(int i=0; i<r_size; i++){double temp1=array[i][i];for(int j=i; j<c_size; j++){array[i][j]/=temp1;}}for(int i=1; i<r_size; i++){for(int j=i-1; j>=0; j--){double temp2=array[j][i];for(int h=i; h<c_size; h++){array[j][h]-=(temp2*array[i][h]);}}}}int main(){int count=0;cin>>count;double **array=new double*[count];for(int i=0; i<count; i++){array[i]=new double[count+1];}for(int i=0; i<count; i++)for(int j=0; j<count+1; j++){cin>>array[i][j];}bool flag=Change(array,count,count+1);if(flag){cout<<endl;double value=Rowlay(array,count,count+1);if(change_num%2==1){value*=-1;}cout<<"行列式为"<<value<<endl;Equation(array,count,count+1);cout<<"方程解为"<<endl;for(int i=0; i<count; i++){char ch='x'+i;cout<<ch<<"= "<<array[i][count]<<endl;}display(array,count,count+1);}else{cout<<"这个矩阵是奇异矩阵"<<endl;cout<<"这个方程无数解！！！"<<endl;cout<<"行列式为"<<0<<endl;}for(int j=0; j<count; j++){delete array[j];}delete []array;}。