数值分析 算法C语言程序
数值分析插值算法源程序
#include<stdio.h>#include<math.h>float f(float x) //计算ex的值{return (exp(x));}float g(float x) //计算根号x的值{return (pow(x,0.5));}void linerity () //线性插值{float px,x;float x0,x1;printf("请输入x0,x1的值\n");scanf("%f,%f",&x0,&x1);printf("请输入x的值: ");scanf("%f",&x);px=(x-x1)/(x0-x1)*f(x0)+(x-x0)/(x1-x0)*f(x1);printf("f(%f)=%f \n",x,px);}void second () //二次插值{float x0,x1,x2,x,px;x0=0;x1=0.5;x2=2;printf("请输入x的值:");scanf("%f",&x);px=((x-x1)*(x-x2))/((x0-x1)*(x0-x2))*f(x0)+((x-x0)*(x-x2))/((x1-x0)*(x1-x2))*f(x1)+((x-x0)* (x-x1))/((x2-x0)*(x2-x1))*f(x2);printf("f(%f)=%f\n",x,px);}void Hermite () //Hermite插值{int i,k,n=2;int flag1=0;printf("Hermite插值多项式H5(x)=");for(i=0;i<=n;i++){int flag=0;flag1++;if(flag1==1){printf("y%d[1-2(x-x%d)*(",i,i);}else{printf("+y%d[1-2(x-x%d)*(",i,i);}for(k=0;k<=n;k++){if(k!=i){flag++;if(flag==1){printf("(1/x%d-x%d)",i,k);}else{printf("+(1/x%d-x%d)",i,k);}}}printf(")]*(");for(k=0;k<=n;k++){if(i!=k){printf("[(x-x%d)/(x%d-x%d)]2",i,k,i);}}printf(")");}printf("\n");}void sectionl () //分段线性插值{float x[5]={2.0,2.1,2.2,2.3,2.4};float y;printf("请输入y:");scanf("%f",&y);if(y>=2.0&&y<2.1){float px;px=((y-x[1])/(x[0]-x[1]))*g (x[0])+((y-x[0])/(x[1]-x[0]))*g (x[1]);printf("f(%f)=%f\n",y,px);}else if(y>=2.1&&y<2.2){float px;px=((y-x[2])/(x[1]-x[2]))*g (x[1])+((y-x[1])/(x[2]-x[1]))*g (x[2]);printf("f(%f)=%f\n",y,px);}else if(y>=2.2&&y<2.3){float px;px=((y-x[3])/(x[2]-x[3]))*g (x[2])+((y-x[2])/(x[3]-x[2]))*g (x[3]);printf("f(%f)=%f\n",y,px);}else if(y>=2.3&&y<2.4){float px;px=((y-x[4])/(x[3]-x[4]))*g (x[3])+((y-x[3])/(x[4]-x[3]))*g (x[4]);printf("f(%f)=%f\n",y,px);}else if(y>2.4) printf("**********ERROR!******************\n"); }void sectionp (){int i;float a[5]={2.0,2.1,2.2,2.3,2.4};float x,y;printf("input the data: x?\n");scanf("%f",&x);if(x<a[1]){i=1;goto loop;}if(x>a[4]){i=4;goto loop;}i=1;loop1:i++;if(x>a[i])goto loop1;if(fabs(x-a[i-1])<=fabs(x-a[i]))i=i-1;loop:y=g(a[i-1])*(x-a[i])*(x-a[i+1])/((a[i-1]-a[i])*(a[i-1]-a[i+1]));y=y+g(a[i])*(x-a[i-1])*(x-a[i+1])/((a[i]-a[i-1])*(a[i]-a[i+1]));y=y+g(a[i+1])*(x-a[i-1])*(x-a[i])/((a[i+1]-a[i-1])*(a[i+1]-a[i]));printf("f(%f)=%f\n",x,y);}int main(){char flag1='y';while(flag1=='y'){int flag=0;printf("*******[1]:线性插值***************\n");printf("*******[2]:二次插值***************\n");printf("*******[3]:Hermite插值************\n");printf("*******[4]:分段线性插值***********\n");printf("*******[5]:分段抛物线插值*********\n");printf("请输入:");scanf("%d",&flag);switch(flag){case 1:linerity ();break;case 2:second ();break;case 3:Hermite ();break;case 4:sectionl ();break;case 5:sectionp ();break;default:printf("error!!\n");}printf("是否继续?y/n \n");getchar();scanf("%c",&flag1);}return 0;}。
数值分析常微分数值解的求法C语言
本科生课程设计报告实习课程数值分析学院名称管理科学学院专业名称信息与计算科学学生姓名学生学号指导教师实验地点实验成绩二〇一六年六月二〇一六年六摘要常微分方程数值解法是计算数学的一个分支.是解常微分方程各类定解问题的数值方法.现有的解析方法只能用于求解一些特殊类型的定解问题,实用上许多很有价值的常微分方程的解不能用初等函数来表示,常常需要求其数值解.所谓数值解,是指在求解区间内一系列离散点处给出真解的近似值.这就促成了数值方法的产生与发展.?关键词:数值解法;常微分1. 实验内容和要求常微分方程初值问题 有精确解2()cos(2)x y x x e x -=+。
要求:分别取步长h = 0.1,0.01,0.001,采用改进的Euler 方法、4阶经典龙格-库塔R -K 方法和4阶Adams 预测-校正方法计算初值问题。
以表格形式列出10个等距节点上的计算值和精确值,并比较他们的计算精确度。
其中多步法需要的初值由经典R-K 法提供。
运算时,取足以表示计算精度的有效位。
2. 算法说明2.1函数及变量说明表1 函数及变量定义1、欧拉方法:1()()(,())i i i i y x y x hf x y x +=+ (i=0,1,2,3,......n-1)(0)y a= (其中a 为初值)2、改进欧拉方法:~1~111()()(,())()()[(,())(,())]2(0)i i i i i i i i i i y x y x hf x y x hy x y x f x y x f x y x y a ++++=+=++=(i=0,1,2......n-1) (其中a 为初值)3、经典K-R 方法: 11213243()6(,)(,)22(,)22(,)i i i i i i i i i i h y y K f x y h hK f x y K h h K f x y K K f x h y hK +⎧=+⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩4、4阶adams 预测-校正方法 预测: 校正:Adsms 内插外插公式联合使用称为Adams 预测-校正系统,利用外插公式计算预测,用内插公式进行校正。
C语言数据分析统计方法和数据可视化
C语言数据分析统计方法和数据可视化随着信息时代的不断发展,大量的数据被产生和积累。
为了从这些海量数据中提取有价值的信息,数据分析和统计方法成为了必不可少的工具。
在编程领域中,C语言作为一种高效且功能强大的语言,提供了各种数据分析和统计的方法,同时结合数据可视化技术,可以帮助程序员更好地理解和解释数据。
一、基本数据分析统计方法C语言提供了许多基础的数据分析和统计方法,如求和、平均值、中位数、极值等。
通过使用循环和条件语句,可以编写简洁高效的代码来实现这些统计功能。
以下是一些常见的基本统计方法的示例代码:1. 求和:```cdouble calculateSum(double array[], int length) {double sum = 0;for (int i = 0; i < length; i++) {sum += array[i];}return sum;}```2. 平均值:```cdouble calculateAverage(double array[], int length) { double sum = calculateSum(array, length);double average = sum / length;return average;}```3. 中位数:```cdouble calculateMedian(double array[], int length) { if (length % 2 == 0) {return (array[length/2 - 1] + array[length/2]) / 2.0; } else {return array[length/2];}}```二、高级数据分析统计方法除了基本的统计方法外,C语言还支持更加高级的数据分析和统计方法,如方差、标准差、相关性等。
这些方法能更全面地描述数据的分布和关系。
1. 方差和标准差:方差度量了数据集合的离散程度,标准差是方差的平方根。
数值计算C语言常用小程序
数值计算C语言常用小程序C语言是一种流行的编程语言,广泛应用于计算机科学和软件开发领域。
在数值计算方面,C语言提供了一些常用的函数和技巧,可以帮助我们进行各种数值计算任务。
本文将介绍一些C语言常用的数值计算小程序。
1.求平均数```c#include <stdio.h>int maiint num;int sum = 0;int count = 0;float average;printf("请输入数字: ");while(num != -1)scanf("%d", &num);if(num != -1)sum += num;count += 1;}}average = (float)sum / count;printf("平均值为: %.2f\n", average);return 0;```这个程序会要求用户输入一系列数字,直到输入-1为止。
然后计算这些数字的平均值并输出。
2.求阶乘```c#include <stdio.h>int factorial(int n)if(n == 0)return 1;} elsereturn n * factorial(n-1);}int maiint n;int result;printf("请输入一个正整数: ");scanf("%d", &n);result = factorial(n);printf("%d的阶乘为: %d\n", n, result);return 0;```这个程序会要求用户输入一个正整数,然后使用递归的方式计算该整数的阶乘,并输出结果。
3.求平方根```c#include <stdio.h>#include <math.h>int maidouble num;printf("请输入一个数字: ");scanf("%lf", &num);if(num < 0)printf("无法计算负数的平方根\n");} elsedouble result = sqrt(num);printf("该数字的平方根为: %.2lf\n", result);}return 0;```这个程序会要求用户输入一个数字,然后计算该数字的平方根并输出。
C++编写数值分析程序
高斯消去法;#include <iostream>using namespace std;template <class T>T** Allocation2D(int m, int n){T **a;a = new T*[m];for (int i=0; i<m; i++){a[i] = new T[n];}return a;}int main(){int i,j,k;int n;float** a;cout<<"输入系数矩阵的N值,N:";cin>>n;a = Allocation2D<float>(n, n+1);cout<<endl<<"输入增广矩阵的各值:/n";for(i=0; i<n; i++){for(j=0; j<n+1; j++){cin>>a[i][j];}}for(k=0; k<n-1; k++){for(i=k+1; i<n; i++){for(j=k+1; j<n+1; j++){a[i][j] = a[i][j] - a[i][k] * a[k][j] / a[k][k];}}}float temp;a[n-1][n] = a[n-1][n] / a[n-1][n-1];for(k=n-2; k>=0; k--){temp = 0;for(j=k+1; j<n; j++){temp = temp + a[k][j] * a[j][n];}a[k][n] = (a[k][n] - temp) / a[k][k];}for(i=0; i<n; i++){cout<<"x"<<i<<": "<<a[i][n]<<endl;}return 0;}拉格朗日插值法#include<iostream>#include<math.h>#define N 3using namespace std;double fun(double *x,double *y, int n,double p);void main(){int i;double a[N],p;double b[N];cout<<"please input xiangliang a= "<<endl;for(i=0;i<N;i++){cin>>a[i];}cout<<"please input xiangliang b= "<<endl;for(i=0;i<N;i++){cin>>b[i];}cout<<"please input LagelangrichazhiJieDian p= "<<endl;cin>>p;cout<<"The Answer= "<<fun(a,b,N,p)<<endl;}double fun(double x[],double y[], int n,double p){double z=0,s=1.0;int k,i=0;double L[N];k=0;while(k<n){if(k==0){for(i=1;i<n;i++){s=s*(p-x[i])/(x[0]-x[i]);}L[0]=s*y[0];k=k+1;}else{s=1.0;for(i=0;i<=k-1;i++){s=s*((p-x[i])/(x[k]-x[i]));}for(i=k+1;i<n;i++){s=s*((p-x[i])/(x[k]-x[i]));}L[k]=s*y[k];k++;}}for(i=0;i<n;i++){z=z+L[i];}return z;}Newton插值法;#include <iostream.h>#include <math.h>void main(){char L;do{double M[100][100];double x[100],y[100];double X=1,xx=0,w=1,N=0,P,R=1;int n;cout<<"请输入所求均差阶数:";cin>>n;for(int i=0;i<=n;i++){cout<<"请输入x"<<i<<"的值:"<<endl;cin>>x[i];cout<<"请输入y"<<i<<"的值:"<<endl;cin>>y[i];M[i][0]=x[i];M[i][1]=y[i];}for( int j=2;j<=n+1;j++){for( i=1;i<=n;i++){M[i][j]=(M[i][j-1]-M[i-1][j-1])/(M[i][0]-M[i-j+1][0]);}}for(i=1;i<=n;i++){cout<<"其"<<i<<"阶均差为:"<<M[i][i+1]<<endl;}cout<<"请输入x的值:x=";cin>>xx;for(i=0;i<n;i++){X*=xx-x[i];N+=M[i+1][i+2]*X;P=M[0][1]+N; cout<<"其函数值:y="<<P<<endl;cout<<endl<<"如还想算其它插值请按'y'否则按'n'"<<endl;cin>>L;}while(L=='y'); }高斯列主元消去法;#include<iostream>#include<cmath>using namespace std;int main(){double max(double a, double b,double c=0);int i,j,t;double temp[3];double a[3][3]={0.012,0.01,0.167,1,0.8334,5.91,3200,1200,4.2}; double b[3]={0.6781,12.1,981};double x[3];for (i=0;i<=2;i++){if(a[i][0]=max(a[0][0] ,a[1][0],a[2][0]))t=i;for(i=0;i<=2;i++)temp[i]=a[t][i];a[t][i]=a[0][i];a[0][i]=temp[i];}double l21=a[1][0]/a[0][0];a[1][1]=a[1][1]-a[0][1]*l21;a[1][2]=a[1][2]-a[0][2]*l21;b[1]=b[1]-b[0]*l21;double l31=a[2][0]/a[0][0];a[2][1]=a[2][1]-a[0][1]*l31;a[2][2]=a[2][2]-a[0][2]*l31;b[2]=b[2]-b[0]*l31;for (i=0;i<=2;i++){if(a[i][1]==max(a[1][0] ,a[2][0]))t=i;for(i=0;i<=2;i++)temp[i]=a[t][i];a[t][i]=a[0][i];a[0][i]=temp[i];}double l32=a[2][1]/a[1][1];a[2][2]=a[2][2]-a[1][2]*l32;b[2]=b[2]-b[1]*l32;x[2]=b[2]/a[2][2];x[1]=(b[1]-a[1][2]*x[2])/a[1][1];x[0]=(b[0]-a[0][1]*x[1]-a[0][2]*x[2])/a[0][0];cout<<"x1="<<x[0]<<" "<<"x2="<<x[1]<<" "<<"x3="<<x[2]<<endl; return 0;}double max(double a,double b,double c){if(a>b){if(b>c) return a;elseif(a>c) return a; else return c; }elseif(b>c) return b; else return c;}迭代法;#include<iostream>using namespace std;double f(double);int main(){double a,i=0,a1,pr;cout<<"please input the start number:"; cin>>a;cout<<"please input precision:";cin>>pr;do{a1=a;a=f(a);i++;if(i>1000000){cout<<"Iterative is defeated!";return 0;} }while((a-a1)>pr||(a-a1)<-pr);cout<<"Iterative Numbers is "<<i<<endl; cout<<"The answer is "<<a<<endl; double f(double x){return pow((x+1.0),(1.0/3));}数值积分(梯形和辛普生公式)#include<stdio.h>double f(double x){return(x/(4+x*x));}double S(double a,double b){double u,h,c;h=(b-a)/6;c=(a+b)/2;u=f(a)+f(b)+4*f(c);return (h*u);}double F(double x){return(x/(4+x*x));}double T(double a,double b){double u,h;h=(b-a)/2;u=F(a)+F(b);return (h*u);}void main(){int c; //循环控制量int i=1; //求解次数printf("----数值积分求解-----\n");printf("1、利用辛普生公式;\n");printf("2、利用梯形公式;\n");printf("0、退出。
数值分析程序
拉格朗日#include<iostream>#include <cmath>using namespace std;double func (double m);int main(){const int j=100;int nn=0;double x=0.0,y=0.0,temp=0;cout<<"请输入要插节点个数"<<endl;cin>>nn;cout<<"请输入x的值:"<<endl;cin>>x;double a[j],b[j];for(int i=0;i<nn;++i){cout<<"请输入x的第"<<i+1<<"个值:";cin>>a[i];temp=a[i];b[i]=func(temp);}for(int k=0;k<nn;++k){double t=1.0;for(int n=0;n<nn;++n){if(k==n) ++n;if(n==nn) break;t=t*(x-a[n])/(a[k]-a[n]);}y=y+t*b[k];}cout<<"结果为:"<<y<<endl;return 0;}double func (double m){return sqrt(m) ;}四阶龙格库塔#include<iostream>using namespace std; double g(double,double);double rungekutta(double x0,double y0,double h,double N);int main(){double x0,y0,h;int N;cout<<"Please input the beginning point x0, y0, step length h, calculating times N."<<endl;cin>>x0>>y0>>h>>N;cout<<rungekutta(x0,y0,h,N);}double rungekutta(double x0,double y0,double h,double N){int i;double x1,y1,K1,K2,K3,K4;for(i=0;i!=N;i++){x1=x0+h;K1=g(x0,y0);K2=g(x0+h/2,y0+h/2*K1);K3=g(x0+h/2,y0+h/2*K2);K4=g(x1,y0+h*K3);y1=y0+h/6*(K1+2*K2+2*K3+K4);x0=x1;y0=y1;}return y1;}//注:以下的函数需要根据实际情况进行调整double g(double x,double y){return y-2*x/y;}快速弦截法#include<iostream>using namespace std;//下列函数的参数表按实际函数的不同定义一般只需要xdouble h(double a,double b,double c,double x);int main(){double x,x1,x2,eps,n,a,b,c;int k=1;//此处求解的方程为a*x^2+b*x+c=0,根据情况改变cout<<"Solving Equation a*x^2+b*x+c=0, input a,b,c,x0,x1,iteration times N and precision Eps.\n";//根据函数本有的变量进行输入一般只需要输入第一个起始点x 第二个起始点x1 次数n和精度epscin>>a>>b>>c>>x>>x1>>n>>eps;while(k<=n){if(h(a,b,c,x1)-h(a,b,c,x)==0){cout<<"Iteration Failed. Try other starting points.\n";return 0;}x2=x1-(x1-x)*h(a,b,c,x1)/(h(a,b,c,x1)-h(a ,b,c,x));if(x2-x1<eps&&x1-x2<eps){cout<<"The solution is x="<<x2<<endl;return 0;}x=x1;x1=x2;++k;}cout<<"Iteration Failed. Cannot be precise enough.\n";return 0;}//以下根据情况改变一般只需要未知数x double h(double a,double b,double c,double x)//Equation ax^2+bx+c=0{return (a*x*x+b*x+c);} 高斯-赛德尔(有点长)#include<iostream>using namespace std;int main(){int N,count,k,n;k=1;cout<<"Please input the number of equations.\n";cin>>N;cout<<"Please input max iretating number.\n";cin>>n;double **a=new double*[N];for(int i=0;i<N;++i){a[i]=new double [N];}double *b=new double[N];double *x=new double[N];double *y=new double[N];//上边的部分最好按自己的习惯定义数组此处给出的是可以随意调整方程阶数的解决方案double eps;cout<<"Please input the factors of x1,x2.....xn.\n";for(int i=0;i<N;++i){for(int j=0;j<N;++j){cin>>a[i][j];}x[i]=0;y[i]=x[i];}cout<<"Please input every result of the equations.\n";for(int i=0;i<N;++i){cin>>b[i];}cout<<"Please input the precision required.\n";cin>>eps;double sum;while(k<=n){for(int i=0;i<N;++i){sum=0.0;for(int j=0;j<N;++j){if(j==i)continue;sum+=a[i][j]*y[j];}y[i]=(b[i]-sum)/a[i][i];}count=0;for(int i=0;i<N;++i){if(x[i]-y[i]<=eps&&y[i]-x[i]<=eps){count++;}}cout<<"k="<<k-1;for(int i=0;i<N;++i){cout<<"y"<<i+1<<'='<<y[i]<<" ";}cout<<endl;if(count==N){cout<<"The result is:\n";for(int i=0;i<N;++i){cout<<"y"<<i+1<<'='<<y[i]<<" ";}return 0;}++k;for(int i=0;i<N;i++)x[i]=y[i];}cout<<"Iteration failed. Can not reach destination precision."<<endl;for(int i=0;i<N;++i){delete[N] a[i];}return 0;}高斯消去法#include<iostream>using namespace std;int main(){double **a;double *b;double *x,*y;int i,j,N;i=1;cout<<"Please input the number of the equations."<<endl;cin>>N;a=new double*[N+1];b=new double[N+1];x=new double[N+1];y=new double[N+1];for(i=1;i<=N;i++){a[i] = new double[N+1];}//以上按习惯定义数组cout<<"Please input the factors of x1,x2.....xn."<<endl;for(i=1;i<=N;i++)for(j=1;j<=N;j++){cin>>a[i][j];}cout<<"Please input the result of every equation."<<endl;for(i=1;i<=N;i++){cin>>b[i];}int k=1;int m;double sum,max,temp;i=1;do{max=fabs(a[k][k]);m=k;for(j=k;j<=N;j++){if(max<fabs(a[j][k])){max=a[j][k];m=j;}if(max==0){cout<<"No solution or Unlimited Solutions. "<<endl;return 0;}}for(int p=1;p<=N;p++){temp=a[m][p];a[m][p]=a[k][p];a[k][p]=temp;}temp=b[m];b[m]=b[k];b[k]=temp;for(j=k+1;j<=N;j++){a[k][j]=a[k][j]/a[k][k];}b[k]=b[k]/a[k][k];for(i=k+1;i<=N;i++){for(j=k+1;j<=N;j++){a[i][j]=a[i][j]-a[i][k]*a[k][j];}b[i]=b[i]-a[i][k]*b[k];}}while(k++!=N);for(i=N-1;i>=1;i--){sum=0;for(j=i+1;j<=N;j++){sum+=a[i][j]*b[j];}b[i]=b[i]-sum;}cout<<"The result is: "<<endl;for(i=1;i<=N;i++)cout<<'x'<<i<<'='<<b[i]<<" ";delete a;delete b;delete x;delete y;return 0;}高斯消去法#include<iostream>#include<cmath>const int n=3;using namespace std;int main(){int k=0;double a[n][n],b[n];cout<<"请输入各个系数:\n";for(int i=0;i<n;++i)for(int j=0;j<n;++j)cin>>a[i][j];cout<<"请输入各方程右边的值\n";for(int j=0;j<n;++j)cin>>b[j];do{double max=fabs(a[k][k]);for(int l=k;l<n;++l)if(max<fabs(a[l][k]) ) max=fabs(a[l][k]);for(int ll=k;ll<n;++ll)if(max==fabs(a[ll][k])){double temp;for(int i=k;i<n;++i){temp=a[k][i];a[k][i]=a[ll][i];a[ll][i]=temp;}double mm=b[k];b[k]=b[ll];b[ll]=mm;}for(int j=k+1;j<n;j++)a[k][j]=a[k][j]/a[k][k];b[k]=b[k]/a[k][k];for(int jj=k+1;jj<n;++jj){int i=jj;a[i][jj]=a[i][jj]-a[i][k]*a[k][jj];}for(int i=k+1;i<n;++i)b[i]=b[i]-a[i][k]*b[k];++k;}while(k!=n-1);for(int ii=n-2;ii>=0;--ii){double sum=0.0;for(int jjj=ii+1;jjj<n;++jjj)sum=sum+a[ii][jjj]*b[jjj];b[ii]=b[ii]-sum;}for(int mn=0;mn<n;++mn)cout<<"x["<<mn<<"]="<<b[mn]<<endl;return 0;}。
数值分析实验报告程序
一、Newton插值法#include<stdio.h>#define MAX_N 20typedef struct tagPOINT{double x;double y;}POINT;int main(){int n,i,j;POINT points[MAX_N+1];double diff[MAX_N+1]; double x,tmp,newton=0;printf("\nInput n value:");scanf("%d",&n);printf("Now input the (x_i,y_i),....,%d:\n",n);for(i=0;i<=n;i++)scanf("%lf%lf",&points[i].x,&points[i].y);printf("Now input the x value:");scanf("%lf",&x);for(i=0;i<=n;i++) diff[i]=points[i].y;for(i=0;i<=n;i++){for(j=n;j>i;j--){diff[j]=(diff[j]-diff[j-1])/(points[j].x-points[j-1-i].x); }}tmp=1;newton=diff[0];for(i=0;i<n;i++){tmp=tmp*(x-points[i].x);newton=newton+tmp*diff[i+1];}printf("newton(%f)=%f\n",x,newton);return 0;}二、Lagrange插值法double Lagrange(int n,double u,double x[],double y[]) {int i,j;double l,Ln=0;for(i=0;i<=n;i++){l=1.0;for(j=0;j<=n;j++)if (i!=j)l=l*(u-x[j])/(x(i)-x[j]);Ln=Ln+l*y[i];}return(Ln);}#include<stido.h>#define N 20V oid main(void){int i,n;double x[N],y[N],u.fu;double Lagrange(int n,double u,double x[],double y[]); printf(“请输入插值多项式的次数N及插入点U:\n”); scanf(“%d,%lf”,&n,&u);printf(“请输入n+1个插值节点x,y:\n”);for(i=0;i<=n;i++)scanf(“%lf,%lf”,&x[i],&y[i]);fu=Lagrange(n,u,x,y);printf(“插入点u的近似值为%lf\n”,fu);}三、用梯形公式#include<stdio.h>#include<math.h>double f(double x){return(sqrt(pow(x,4)+1));}double T(double a,double b){double u,h;h=(b-a)/2;u=f(a)+f(b);return (h*u);}void main(){ double a,b;printf("请输入积分区域:a,b\n"); scanf("%lf,%lf",&a,&b);double T(double ,double );printf("%lf,%lf",T(a,b));四、Simpson公式求积分#include<stdio.h>#include<math.h>double f(double x){return(sqrt(2-pow(sin(x),2))); }double S(double a,double b){double u,h,c;h=(b-a)/6;c=(a+b)/2;u=f(a)+f(b)+4*f(c);return (h*u);}void main(){ double a,b;printf("请输入积分区域:a,b\n");scanf("%lf,%lf",&a,&b);double S(double ,double );printf("%lf ",S(a,b))五、复化数值积分公式的自动控制误差算法复化Simpson自动控制误差#include<math.h>double AS(double (*f)(double x),double a,double b,double e) {int i;int n=2;double h,Jn,J2n,Sn,S2n;h=(b-a)/n;Jn=(*f)((a+b)/2);Sn=h/3*((*f)(a)+(*f)(b)+4*Jn);while(1){J2n=0;for(i=1;i<=n;i++)J2n=J2n+(*f)(a+(2*i-1)*h/2);S2n=Sn/2+h/6*(4*J2n-2*Jn);if(fabs(S2n-Sn)<15*e){return S2n;break;}else{ h=h/2;n=2*n;J2n=Jn;Sn=S2n;}}}复化梯形公式double comTx(double (*f)(double),double a,double b,int n) {int i; double h,Tn;h=(b-a)/n;Tn=(*f)(a)+(*f)(b);for(i=1;i<=n-1;i++)Tn=Tn+2*(*f)(a+i*h);Tn=(h/2)*Tn;return(Tn);}复化梯形自动控制误差double A T(double (*f)(double x),double a,double b,double e) {int i; int n=10;double h,Js,Tn,T2n;h=(b-a)/n;Tn=comTx(f,a,b,n);while(1){Js=0;for(i=1;i<=n;i++)Js=Js+(*f)(a+(2*i-1)*h/2);T2n=Tn/2+h/2*Js;if(fabs(T2n-Tn)<3*e)break;else{ h=h/2;n=2*n;Tn=T2n;}}return T2n;}六、Romberg公式计算积分#include "stdio.h"#include "math.h"#define MAX 20double f(double x)//要求的积分公式{ return(lnx);}double computeTn(double (*f)(double x),double a,double b, int n)//复化梯形公式{int i;double sum=0;double h=(b-a)/n;for(i=1;i<n;i++)sum+=f(a+i*h);sum+=(f(a)+f(b))/2;return (h*sum);}double Romberg(double (*f)(double x),double a,double b,double ep)//Romberg求积分{ int j,n=1,k;double R[3][MAX];R[2][1]=computeTn(f,a,b,n);printf("%lf\n",R[2][1]);for(k=2;k<MAX;k++){ for(j=1;j<k;j++)R[1][j]=R[2][j];n*=2;R[2][1]=computeTn(f,a,b,n);printf("%lf\t",R[2][1]);for(j=2;j<=k;j++){ R[2][j]=R[2][j-1]+(R[2][j-1]-R[1][j-1])/(pow(4,j-1)-1);printf("%lf\t",R[2][j]);}printf("\n");if(fabs(R[1][k-1]-R[2][k])<ep) return(R[2][k]); }k--;if(fabs(R[1][k-1]-R[2][k])>ep){printf("Return no solved...\n");return 0;}}七、用牛顿迭代法、弦截法求非线性方程的根///用牛顿迭代法的C源程序///#include<stdio.h>#include<math.h>#define m 20double f(double x){return((x*(x*(x-7.7)+19.2)-15.3));}double f1(double x){return(x*(3*x-15.4)+19.2);}main(){int k;double x,x0,e;printf("请输入x0和允许误差e的值:\n");scanf("%lf,%lf",&x0,&e);for(k=1;k<=m;k++){x=x0-f(x0)/f1(x0);if(fabs(x-x0)<e){printf("方程的近似根为%lf",x);return NULL;} x0=x;}printf("给定的迭代次数太小.停机");}///用弦截法的C源程序///#include<stdio.h>#include<math.h>#define m 20double f(double x){return((x*(x*(x-7.7)+19.2)-15.3));}main(){int k;double x,x0,x1,f0,f1,e;printf("请输入x0,x1和允许误差e的值:\n");scanf("%lf,%lf,%lf",&x0,&x1,&e);f0=f(x0);for(k=1;k<=m;k++){f1=f(x1);x=x1-(f1*(x1-x0)/(f1-f0));if(fabs(x-x1)<e){printf("方程的近似根为%lf",x);return NULL;} f0=f1;x0=x1;x1=x;}printf("给定的迭代次数太小.停机");}八、Gauss-Seidel迭代法解线性方程组的算法#include <math.h>#include <stdio.h>#define eps 0.000001#define max 100#define N 20double norm_inf(double x[],int n){ double norm;int i;norm=fabs(x[0]);for(i=1;i<n;i++){if(fabs(x[i])>norm)norm=fabs(x[i]);}return norm;}void seidel(double a[N][N],double g[N],int n){double b[N][N]={0},x0[N],x1[N],x1_x0[N],norm,temp; int i,j,k;for(i=0;i<n;i++){g[i]=g[i]/a[i][i];for(j=0;j<n;j++){ if(i==j)continue;b[i][j]=-a[i][j]/a[i][i];}}for(i=0;i<n;i++){x0[i]=0;x1[i]=1;x1_x0[i]=x1[i]-x0[i];}k=0;norm=norm_inf(x1_x0,n); while((norm>=eps)&&(k<max)) {for(i=0;i<n;i++)x0[i]=x1[i];for(i=0;i<n;i++){temp=0;for(j=0;j<=i-1;j++)temp=temp+b[i][j]*x1[j];for(j=i+1;j<n;j++)temp=temp+b[i][j]*x0[j];x1[i]=temp+g[i];x1_x0[i]=x1[i]-x0[i];}norm=norm_inf(x1_x0,n);k++;}for(i=0;i<n;i++)printf("x[%d]=%lf\n",i,x1[i]); printf("%d times iteration.\n",k); }int main(){double b[N][N]={{10,-2,-1},{-2,10,-1},{-1,-2,5}},f[N]={0,-21,-20}; printf("\n");seidel(b,f,3);return 0;}九、用列主元消元法#include<stdio.h>#include<math.h>#define n 3void print(double a[n][n+1]);void gauss(double a[n][n+1],double x[n]){ int i,j,k;double temp,s,l;for(i=0;i<n-1;i++){ k=i;for(j=i+1;j<n;j++){ if(fabs(a[j][i])>fabs(a[k][i]))k=j;}if(k!=i)for(j=i;j<=n;j++) {temp=a[i][j];a[i][j]=a[k][j];a[k][j]=temp;}for(j=i+1;j<n;j++){ l=1.0*a[j][i]/a[i][i];for(k=0;k<n+1;k++)a[j][k]=a[j][k]-a[i][k]*l; }print(a);printf("lf");}print(a);x[n-1]=a[n-1][n]/a[n-1][n-1]; for(i=n-2;i>=0;i--){ s=0.0;for(j=i;j<n;j++){ if(j==i)continue;s+=a[i][j]*x[j];}x[i]=(a[i][n]-s)/a[i][i];}}void print(double a[n][n+1]){int i,j;for(i=0;i<n;i++){ for(j=0;j<n+1;j++)printf("%lf",a[i][j]);printf("lf");}}十、Doolittle分解法#include "iostream.h"#include "stdio.h"#define max 20void main(){int i,k,n,j,p;double x,v,y;double a[max][max],r[max][max],l[max][max]; double b[max],Y[max],X[max];printf("Input n:"); scanf("%d",&n);printf("Input a[%d][%d]:\n",n,n);for(i=1;i<n+1;i++)for(j=1;j<n+1;j++){ scanf("%lf",&x); a[i][j]=x; }printf("b[i]:");for(i=1;i<n+1;i++){ scanf("%lf",&y); b[i]=y; }for(k=1;k<n+1;k++){ for(i=k;i<n+1;i++){ v=0;for(p=1;p<k;p++)v=l[k][p]*r[p][i]+v;r[k][i]=a[k][i]-v; printf("r[%d][%d]=%f\n",k,i,r[k][i]); } for(i=k+1;i<n+1;i++){ v=0;for(p=1;p<k;p++)v=v+l[i][p]*r[p][k];l[i][k]=(a[i][k]-v)/r[k][k]; printf("l[%d][%d]=%f\n",i,k,l[i][k]);} for(i=1;i<n+1;i++){ v=0;for(k=1;k<i;k++)v=v+l[i][k]*Y[k];Y[i]=b[i]-v;}for(i=n;i>0;i--){ v=0;for(k=i+1;k<n+1;k++)v=v+r[i][k]*X[k];X[i]=(Y[i]-v)/r[i][i];}printf("\n");for(i=1;i<n+1;i++)printf("x[%d]=%f\n",i,X[i]);printf("\n");return ;}。
数值分析C程序
b[k*n+i] = t;
}
}
// 为了使a的第k行的第k元素为1,a和b的第k行都除以a(k,k)
t = 1.0 / a[k*n+k];
for(i=k; i<n; i++) a[k*n+i] *= t;
for(i=0; i<n; i++) b[k*n+i] *= t;
for(i=0; i<n; i++)
{
if(i != k && a[i*n+k] != 0.0)
{
t = a[i*n+k];
for(j=k+1; j<n; j++) a[i*n+j] -= a[k*n+j] * t;
for(j=0; j<n; j++) b[i*n+j] -= b[k*n+j] * t;
// 带有列主元的约当消元法
// 功能: 求解线性方程组 Ax = b
// 参数: A - 指向n*n系数矩阵的指针
// b - 常数向量的指针
// n - 方程组的维数
// 返回值:0 - 如果成功。线性方程组的解保存在 b 中
// 1 - 求解失败
// ----------------------------------------------------------------------------
if(j != k) // 交换方程的第 j 行和 k 行
{
for(i=k; i<n; i++)
数值分析五个题目的C语言及Matlab程序
本文档包含上一个文档中的五个数值分析实验题C语言程序及Matlab程序实验一C程序#include "stdio.h"#include "math.h"void main(){inti=0;float a=0.1,b=1.9,t=0.0,e=1.9;if((pow(a,7)-28*pow(a,4)+14)*(pow(b,7)-28*pow(b,4)+14)<0) if((7*pow(x,6)-112*pow(x,3)))printf("x=%f,i=%d,e=%f\n",x,i,e);for(i=1;i<7&&e>0.00001;i++){t=x;x=x-(pow(x,7)-28*pow(x,4)+14)/(7*pow(x,6)-112*pow(x,3));e=fabs(t-x);printf("x=%f,i=%d,e=%f\n",x,i,e);}}Matable 程序i=0;x=1.9;t=0.0;e=1.9;disp(['i=',num2str(i),' ','x=',num2str(x),' ','e=',num2str(e)]); for i=1:7t=x;x=x-(x^7-28*x^4+14)/(7*x^6-112*x^3);e=abs(t-x);disp(['i=',num2str(i),' ','x=',num2str(x),' ','e=',num2str(e)]);if e<0.00001break;endend实验二C程序#include"stdio.h"#include"math.h"//已知量double x[10]={1,2,3,4,5,6,7,8,9,10};double fd1=1,fd2=0.1;double fx[10]={0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.9459101,2.079445,2.1972246,2.3025851}; //函数声明void canshu(double *h,double *r,double *u,double *d);void wanju(double *u,double *r,double *d,double *m);double yths(double *m,double *h,double p);double ytdhs(double *m,double *h,double p);//主函数void main(){double p,q1,q2;double h[9],u[9],r[9],m[10],d[10];// 矩阵Am=d参数数据u[8]=1,r[0]=1;canshu(h,r,u,d);wanju(u,r,d,m);printf("请输入待求数值:x=");scanf("%lf",&p);q1=yths(m,h,p);q2=ytdhs(m,h,p);printf("输出结果为:\nf(%lf)=%lf\nf'(%lf)=%lf\n",p,q1,p,q2);}//求参函数void canshu(double *h,double *r,double *u,double *d) {inti,j;double a[10];for(i=0;i<9;i++){h[i]=*(x+i+1)-*(x+i);a[i]=(fx[i]-fx[i+1])/(x[i]-x[i+1]);j=i-1;if(j>-1){r[i]=h[j+1]/(h[j+1]+h[j]);u[j]=h[j]/(h[j]+h[j+1]);d[i]=(6/(h[j]+h[i]))*(a[i]-a[j]);}}d[0]=(6/h[0])*(a[i-9]-fd1);d[9]=(6/h[i-1])*(fd2-a[i-1]);}//追赶法求弯矩向量mvoid wanju(double *u,double *r,double *d,double *m){inti,j;double l[9],n[10],z[10];n[0]=2;z[0]=d[0];for(i=0;i<9;i++){l[i]=u[i]/n[i];n[i+1]=2-l[i]*r[i];z[i+1]=d[i+1]-l[i]*z[i];}m[9]=z[9]/n[9];for(j=8;j>-1;j--){m[j]=(z[j]-r[j]*m[j+1])/n[j];}}//三次样条插值函数double yths(double *m,double *h,double p) {inti,j=0;double q;for(i=0;i<10;i++){if(x[i]<=(p-1)&&x[i+1]>=(p-1))break;}q=pow(x[i+1]-p,3)/6/h[i+1]*m[i]+pow(p-x[i],3)/6/h[i+1]*m[i+1]+(fx[i]-pow(h[i+1],2)*m[i]/6)*(x[i+1]-p)/h[i+1]+(fx[i+1]-pow(h[i+1],2)*m[i+1]/6)*(p-x[i])/h[i+1];return q;}//三次样条插值导函数double ytdhs(double *m,double *h,double p){inti,j=0;double q;for(i=0;i<10;i++){if(x[i]<=(p-1)&&x[i+1]>=(p-1))break;}q=-pow(x[i+1]-p,2)/2/h[i+1]*m[i]+pow(p-x[i],2)/2/h[i+1]*m[i+1]- (fx[i]-pow(h[i+1],2)*m[i]/6)/h[i+1]+(fx[i+1]-pow(h[i+1],2)*m[i+1]/6)/h[i+1];return q;}Matlab程序function[h,r,u,d]=canshu(x,fx,h,r,u,d,fd1,fd2)%UNTITLED10 Summary of this function goes here % Detailed explanation goes herefor i=1:9h(i)=x(i+1)-x(i);a(i)=(fx(i)-fx(i+1))/(x(i)-x(i+1));j=i-1;if j>0r(i)=h(j+1)/(h(j+1)+h(j));u(j)=h(j)/(h(j)+h(j+1));d(i)=(6/(h(j)+h(i)))*(a(i)-a(j));endendd(1)=(6/h(1))*(a(i-8)-fd1);d(10)=(6/h(i))*(fd2-a(i));endglobal x fd1 fd2 fx h d m r u;x=[1,2,3,4,5,6,7,8,9,10];fd1=1;fd2=0.1;fx=[0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.9459101,2.079445,2.19 72246,2.3025851];u(9)=1;r(1)=1;[h,r,u,d]=canshu(x,fx,h,r,u,d,fd1,fd2);[m]=wanju(u,r,d,m);p=input('请输入待求数值:x=');q1=yths(m,h,p,x,fx);q2=ytdhs(m,h,p,fx,x);disp(['输出结果为:','f(',num2str(p),')=',num2str(q1),',','f(',num2str(p),')=',num2str(q2)]); function[m]=wanju(u,r,d,m)%UNTITLED12 Summary of this function goes here% Detailed explanation goes heren(1)=2;z(1)=d(1);for i=1:9l(i)=u(i)/n(i);n(i+1)=2-l(i)*r(i);z(i+1)=d(i+1)-l(i)*z(i);endm(10)=z(10)/n(10);for j=9:-1:1m(j)=(z(j)-r(j)*m(j+1))/n(j);endendfunction [q] =ytdhs(m,h,p,fx,x)%UNTITLED14 Summary of this function goes here% Detailed explanation goes herefor i=1:10if x(i)<=(p)&&x(i+1)>=(p)break;endq=-(x(i+1)-p)^2/2/h(i+1)*m(i)+(p-x(i))^2/2/h(i+1)*m(i+1)-(fx(i)-(h(i+1))^2*m(i)/6)/h(i+1)+(fx(i +1)-(h(i+1))^2*m(i+1)/6)/h(i+1);endendfunction[q]=yths(m,h,p,x,fx)%UNTITLED13 Summary of this function goes here% Detailed explanation goes herefor i=1:10if x(i)<=(p)&&x(i+1)>=(p)break;endq=(x(i+1)-p)^3/6/h(i+1)*m(i)+(p-x(i))^3/6/h(i+1)*m(i+1)+(fx(i)-(h(i+1))^2*m(i)/6)*(x(i+1)-p)/h(i+1)+(fx(i+1)-(h(i+1))^2*m(i+1)/6)*(p-x(i))/h(i+1); endend实验三C程序#include"stdio.h"#include"math.h"//已知量double b=3,a=1;//函数声明double romberg();double fhs(double x);//主函数void main(){ double q;q=romberg();printf("输出结果为:r=%lf\n",q);}//Romberg算法double romberg(){intn,j,k=0,m=0,z=0,i;double e=1,qh=0,q,h,x0,r1=0,t[50],s[50],c[50],r[50]; t[0]=(b-a)/2*(fhs(a)+fhs(b));for(i=1;1;i++){n=pow(2,i-1);h=(b-a)/n;x0=a+0.5*h;for(j=0;j<n;j++){qh+=fhs(x0+j*h);}t[i]=0.5*(t[i-1]+h*qh);qh=0;if(i>=3){for(k;k<i;k++)s[k]=(4*t[k+1]-t[k])/3;for(m;m<k-1;m++)c[m]=(16*s[m+1]-s[m])/15;for(z;z<m-1;z++)r[z]=(64*c[z+1]-c[z])/63;e=fabs(r1-r[z-1]);r1=r[z-1];}if(e<0.00001)break;}q=r[z-1];return q;}//f函数double fhs(double x){double y;y=pow(3,x)*pow(x,14)*(5*x+7)*sin(pow(x,2));return y;}Matlab程序function [y] =fhs(x)%UNTITLED4 Summary of this function goes here % Detailed explanation goes herey=3^x*x^14*(5*x+7)*sin(x^2);endfunction[q]=romberg()%UNTITLED3 Summary of this function goes here % Detailed explanation goes hereglobal b a n j k m z i e qh h x0;t=zeros(1,20);s=zeros(1,20);c=zeros(1,20);r=zeros(1,20);k=1;m=1;z=1;e=1;r1=0;t(1)=(b-a)/2*(fhs(a)+fhs(b));qh=0;for i=1:20n=2^(i-1);h=(b-a)/n;x0=a+0.5*h;for j=0:n-1qh=qh+fhs(x0+j*h);endt(i+1)=0.5*(t(i)+h*qh);qh=0;if i>=4for k=k:is(k)=(4*t(k+1)-t(k))/3;endfor m=m:k-1c(m)=(16*s(m+1)-s(m))/15;endfor z=z:m-1r(z)=(64*c(z+1)-c(z))/63;ende=abs(r1-r(z-1));r1=r(z-1);endif(e<0.00001)break;endendq=r(z-1);endglobal b a k m z e qh r1;k=0;m=0;z=0;e=1;qh=0;r1=0;b=3;a=1;[q]=romberg();disp(['输出结果为:r=',num2str(q)]);实验四C程序#include"stdio.h"#include"math.h"//已知量double x=0,h=0.0005,y[3]={0},p[4]={0.025,0.045,0.085,0.1}; double fx(double x,double y[3],int t);void main(){inti,t;double k[4],m[3];for(i=0;i<4;i++){for(x=0;x<=p[i];x=x+h){for(t=0;t<3;t++)m[t]=y[t];for(t=0;t<3;t++){k[0]=fx(x,m,t);m[t]=y[t]+0.5*h*k[0];k[1]=fx(x+0.5*h,m,t);m[t]=y[t]+0.5*h*k[1];k[2]=fx(x+0.5*h,m,t);m[t]=y[t]+h*k[2];k[3]=fx(x+h,m,t);//m[t]=y[t];y[t]=y[t]+h*(k[0]+2*k[1]+2*k[2]+k[3])/6;}}x=0;for(t=0;t<3;t++)printf("y%d(%lf)=%lf\t",t+1,p[i],y[t]);printf("\n");for(t=0;t<3;t++)y[t]=0;}}double fx(double x,double y[3],int t){double f;if(t==0)f=0*x+1;else if(t==1)f=0*x+y[2];elsef=0*x+1000-1000*y[1]-100*y[2];return f;}Matlab程序function[f]=fx(x,y)%UNTITLED7 Summary of this function goes here % Detailed explanation goes heref(1)=0*x+1;f(2)=0*x+y(3);f(3)=0*x+1000-1000*y(2)-100*y(3);endy=zeros(1,3);h=0.0005;x=0;p=[0.025,0.045,0.085,0.1];k=zeros(4,3);for i=1:4for j=1:fix(p(i)/h)k(1,:)=fx(x,y);k(2,:)=fx(x+0.5*h,y+0.5*h*k(1,:));k(3,:)=fx(x+0.5*h,y+0.5*h*k(2,:));k(4,:)=fx(x+h,y+h*k(3,:));y=y+h*(k(1,:)+2*k(2,:)+2*k(3,:)+k(4,:))/6;x=x+h;enddisp(['x=',num2str(p(i)),' ','y=',num2str(y)]);y=zeros(1,3);end实验五C程序#include"stdio.h"#include"math.h"#define N 9//已知量doubleA[N][N]={12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.0678 13,-2.031743,2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124,-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103,1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0 .037585,-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1 .784317,0.718719,1.563849,1.836742,-3.741856,0.431637,9.789356,-0.103458,-1.103456,0.238 417,1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.7138465,3.123789,-2. 213474,3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782,-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00 001};doubleb[9]={2.1874369,33.992318,-25.173417,-0.84671695,1.784317,-86.612343,1.1101230,4. 719345,-5.6784392};//函数声明void change(inta,int b);//主函数void main(){ inti,j,q1,n,m,t;double tmp1,tmp2=0,x[9],tmp3;for(i=0;i<9;i++){tmp1=A[i][i];q1=i;for(j=i;j<9;j++)//{if(fabs(tmp1)<fabs(A[j][i])){ tmp1=A[j][i];q1=j;}}if(q1!=i)change(i,q1);for(n=i+1;n<9;n++){tmp3=A[i][i]/A[n][i];for(m=i+1;m<9;m++)A[n][m]=A[n][m]*tmp3-A[i][m];//需检查,为何主元用公式置零不行。
数值分析五个题目的C语言及MATLAB程序
1.7917595,1.9459101,2.079445,2.1972246,2.3025851}; //函数声明 void canshu(double *h,double *r,double *u,double *d); void wanju(double *u,double *r,double *d,double *m); double yths(double *m,double *h,double p); double ytdhs(double *m,double *h,double p); //主函数 void main() {
break; end end 实验二 C 程序 #include"stdio.h" #include"math.h" //已知量 double x[10]={1,2,3,4,5,6,7,8,9,10}; double fd1=1,fd2=0.1; double fx[10]={0,0.69314718,1.0986123,1.3862944,1.6094378,
} Matlab 程序 function[h,r,u,d]=canshu(x,fx,h,r,u,d,fd1,fd2) %UNTITLED10 Summary of this function goes here % Detailed explanation goes here for i=1:9
h(i)=x(i+1)-x(i); a(i)=(fx(i)-fx(i+1))/(x(i)-x(i+1)); j=i-1;
end end function [q] =ytdhs(m,h,p,fx,x) %UNTITLED14 Summary of this function goes here % Detailed explanation goes here for i=1:10
(完整word版)数值分析作业(C语言编程实现)(word文档良心出品)
#include <stdio.h>#include <math.h>double f(double x){double ans;ans=exp(x);return ans;}void main(){double a=1,b=3,error=0.0001,t[20][20],h,c;int i,j,k,m,n;h=b-a;t[0][0]=h*(f(a)+f(b))/2;k=1;while(1){t[0][k]=0;m=1;for(j=0;j<k-1;j++)m=m*2;for(i=1;i<=m;i++)t[0][k]=t[0][k]+h*f(a+(i-0.5)*h);t[0][k]=(t[0][k]+t[0][k-1])/2;for(j=1;j<=k;j++){ c=1;for(n=0;n<j;n++)c=c*4;t[j][k-j]=(c*t[j-1][k-j+1]-t[j-1][k-j])/(c-1);}if(fabs(t[k][0]-t[k-1][0])<error){ printf("\n积分结果I ≈%lf\n",t[k][0]);break;}else{ h=h/2;k++;}}}#include <stdio.h>#include <math.h>double f(double t){double ans;ans=pow(cos(t),1.0/3);return ans;}void main(){double x=0,eslong=0.000001,x0;int N=20,i;printf("\n近似初值x0 = %lf\n",x);for(i=0;i<N;i++){x0=x;x=f(x);printf(" x%d = %lf\n",i+1,x);if(fabs(x-x0)<eslong)break;}if(fabs(x-x0)<eslong)printf("得到近似结果为x ≈%lf\n\n",x,i);elseprintf("迭代失败\n");}#include <stdio.h>#include <math.h>double a=0,b=1,x,y=0,h=0.1,k1,k2,k3,k4;int i,N;double f(double t,double s){double ans;ans=1+t*t;return ans;}void main(){N=(b-a)/h;x=a;printf("\n 初值为(x0,y0) = ( %.8f , %.8f )\n",x,y);for(i=0;i<N;i++){k1=f(x,y);k2=f(x+h/2,y+h*k1/2);k3=f(x+h/2,y+h*k2/2);k4=f(x+h,y+h*k3);y=y+h*(k1+2*(k2+k3)+k4)/6;x=x+h;printf(" 第%d次输出结果为(x%d,y%d) = ( %.8f , %.8f )\n",i+1,i+1,i+1,x,y);}}#include <stdio.h>void main(){double datax[4]={1.2,2.9,4.6,5.8},datay[10]={14.84,33.71,58.36,79.24},l[3],x=1.5,y;int i,j;y=0;for(i=0;i<=3;i++){l[i]=1;for(j=0;j<i;j++)l[i]=(x-datax[j])/(datax[i]-datax[j])*l[i];for(j=i+1;j<=3;j++)l[i]=(x-datax[j])/(datax[i]-datax[j])*l[i];y=y+datay[i]*l[i];}printf("\n f(x)在x = %f 处的近似值为: y = %f\n",x,y);}#include <stdio.h>void main(){double datay[9]={11.7,14.87,21.44,31.39,44.73,61.46,81.57,105.11,131.91};int m=2,i,j,k;double p,data[9][4],a[3][4],datax[9]={1.2,2.3,3.4,4.5,5.6,6.7,7.8,8.9,10.0};for(i=0;i<9;i++)for(j=1;j<2*m+1;j++){data[i][j]=1;for(k=0;k<j;k++)data[i][j]=data[i][j]*datax[i];}for(i=0;i<m+1;i++){for(j=0;j<m+1;j++){a[i][j]=0;for(k=0;k<9;k++)a[i][j]=a[i][j]+data[k][i+j];}}a[0][0]=9;a[0][m+1]=0;for(i=0;i<9;i++)a[0][m+1]=a[0][m+1]+datay[i];for(i=1;i<m+1;i++){ a[i][m+1]=0;for(j=0;j<9;j++){p=datay[j];for(k=0;k<i;k++)p=p*datax[j];a[i][m+1]=a[i][m+1]+p;}} //生成m+1行,m+2列增广矩阵// for(i=0;i<m+1;i++) //显示方程组//for(j=0;j<m+2;j++){if(j!=m+1){ printf("(%f)a%d ",a[i][j],j);if(j!=m)printf("+ ");}elseprintf("= %f \n",a[i][j]);}for(i=0;i<m;i++) //高斯消去法//{if(a[i][i]!=0){ for(j=i+1;j<m+1;j++){ a[j][i]=a[j][i]/a[i][i];for(k=i+1;k<m+2;k++)a[j][k]=a[j][k]-a[i][k]*a[j][i];}}elsebreak;}if(a[m][m]!=0&&i==m){ a[m][m+1]=a[m][m+1]/a[m][m];for(i=2;i<=m+1;i++){ for(j=1;j<i;j++)a[m+1-i][m+1]=a[m+1-i][m+1]-a[m+1-i][m+1-j]*a[m+1-j][m+1];a[m+1-i][m+1]=a[m+1-i][m+1]/a[m+1-i][m+1-i];}printf("方程组的解为:\n");for(j=0;j<m+1;j++)printf("a%d = %f\n",j,a[j][m+1]);printf("拟合多项式为:\n");printf("P%d(x) = (%f) + (%f)x + (%f)x^2\n",m,a[0][m+1],a[1][m+1],a[2][m+1]);}elseprintf("数据有误!\n");}}列主元素法#include <stdio.h>#include <math.h>void main(){double a[3][4]={1,-2,-1,3,-2,10,-3,15,-1,-2,5,10},mov,comp;int i,j,k,nrow;for(i=0;i<2;i++){comp=fabs(a[i][i]);for(k=i;k<3;k++) //比较绝对值大小并进行主元列交换// if(fabs(a[k][i])>=comp){nrow=k;comp=fabs(a[k][i]);}for(j=0;j<=3;j++){mov=a[i][j];a[i][j]=a[nrow][j];a[nrow][j]=mov;}printf("方程第%d行互换位置后如下\n",i+1);for(j=0;j<3;j++)printf("(%f)x1 + (%f)x2 + (%f)x3 = %f\n",a[j][0],a[j][1],a[j][2],a[j][3]);if(a[i][i]!=0){for(j=i+1;j<3;j++){a[j][i]=a[j][i]/a[i][i];for(k=i+1;k<=3;k++)a[j][k]=a[j][k]-a[i][k]*a[j][i];a[j][i]=0;}printf("方程经%d次消元如下\n",i+1);for(j=0;j<3;j++)printf("(%f)x1 + (%f)x2 + (%f)x3 = %f\n",a[j][0],a[j][1],a[j][2],a[j][3]);}elsebreak;}if(a[2][2]!=0&&i==2){printf("方程化简得\n");for(i=0;i<3;i++)printf("(%f)x1 + (%f)x2 + (%f)x3 = %f\n",a[i][0],a[i][1],a[i][2],a[i][3]);a[2][3]=a[2][3]/a[2][2];for(i=2;i<=3;i++){for(j=1;j<i;j++)a[3-i][3]=a[3-i][3]-a[3-i][3-j]*a[3-j][3];a[3-i][3]=a[3-i][3]/a[3-i][3-i];}printf("方程组的解为:\n");for(j=0;j<3;j++)printf("x%d = %f\n",j+1,a[j][3]);}elseprintf("数据有误!\n");}Jacobi迭代法#include <stdio.h>#include <math.h>void main(){double a[3][7]={{1,-2,-1,3},{-2,10,-3,15},{-1,-2,5,10}},error=0.000001,norm;int N=423,i,j,k;a[0][4]=0,a[1][4]=0,a[2][4]=0;for(i=0;i<3;i++) //把a矩阵转化为b矩阵//{a[i][6]=a[i][i];for(j=0;j<3;j++){a[i][j]=-a[i][j]/a[i][6];}a[i][3]=a[i][3]/a[i][6];a[i][i]=0;}printf("化为b矩阵如下\n");for(i=0;i<3;i++){printf("%f %f %f %f\n",a[i][0],a[i][1],a[i][2],a[i][3]);}for(i=1;i<N;i++){for(j=0;j<3;j++){a[j][5]=0;for(k=0;k<3;k++){a[j][5]=a[k][4]*a[j][k]+a[j][5];}a[j][5]=a[j][5]+a[j][3];}norm=0;for(k=0;k<3;k++)norm=norm+fabs(a[k][4]-a[k][5]);if(norm<error)break;elsefor(k=0;k<3;k++)a[k][4]=a[k][5];}if(norm<error){printf("计算结果为\n");for(i=0;i<3;i++){printf(" x%d = %.3f\n",i+1,a[i][5]);}}elseprintf("迭代失败\n");}题目1#include "stdio.h"#include "math.h"double f(double x){double ans;ans=exp(x);return(ans);}void main(){double a=-1,b=1,error=0.0001,m=1,h,T0,T,F;int k;h=(b-a)/2;T0=h*(f(a)+f(b));while(1){F=0;for(k=1;k<=pow(2.0,m-1);k++)F=F+f(a+(2*k-1)*h);T=T0/2+h*F;if(fabs(T-T0)<error)break;m++;h=h/2;T0=T;}printf("积分结果为I ≈%f\n",T);}#include "stdio.h"double f(double t,double s){double ans;ans=1+t*t;return(ans);}void main(){double a=0,b=1,h=0.2,x0=0,y0=0,x,k1,k2,k3,y;int N,n;N=(b-a)/h;for(n=1;n<=N;n++){x=x0+h;k1=f(x0,y0);k2=f(x0+h/2,y0+h/2*k1);k3=f(x0+h,y0-h*k1+2*h*k2);y=y0+h/6*(k1+4*k2+k3);printf("第%d次输出结果为(%.8f,%.8f)\n",n,x,y);x0=x;y0=y;}}。
数值分析 c语言实现高斯赛德尔解法
数值分析程序设计学院:计算机学院姓名:袁薪洋1.实验目的:1熟练掌握C语言程序设计,编程求解问题。
2.运用高斯-赛德尔迭代公式求解线性方程组。
2.实验内容:用高斯-赛德尔迭代公式求解方程组。
10x1-x2-2x3=7.2-x1+10x2-2x3=8.3-x1-x2+5x3=4.2程序的核心代码:#include"math.h"#include<stdio.h>#define NUMBER 20float A[NUMBER][NUMBER+1] ; float ark;int flag,n;void exchange(int r,int k);float max(int k);void main(){float x[NUMBER]; /*此数组用于存放方程解*/int r,k,i,j;printf("***********************************");printf("\n\n用高斯-赛德尔迭代法解线性方程组\n\n");printf("***********************************");printf("\n\n 请输入方程组的维数:n=");scanf("%d",&n);printf(" \n\n请输入系数矩阵A和向量b:");for(i=1;i<=n;i++){printf("\n\n请输入a%d1--a%d%d系数和向量b%d(数之间用空格格开):",i,i,n,i);//实现将每一行中的系数和向量一次性输入,数之间用空格格开,输完后回车确定for(j=1;j<=n+1;j++) //将刚才输入的数存入数组scanf("%f",&A[i][j]);}for(k=1;k<=n-1;k++){ark=max(k);if(ark==0) //判断方程是否为线性方程printf("\n\n此方程组不合法!");}else if(flag!=k)exchange(flag,k);for(i=k+1;i<=n;i++)for(j=k+1;j<=n+1;j++)A[i][j]=A[i][j]-A[k][j]*A[i][k]/A[k][k]; }x[n]=A[n][n+1]/A[n][n];for( k=n-1;k>=1;k--){float me=0;for(j=k+1;j<=n;j++){me=me+A[k][j]*x[j];}x[k]=(A[k][n+1]-me)/A[k][k];}for(i=1;i<=n;i++){printf(" \n\nx%d=%f",i,x[i]);}void exchange(int r,int k) //交换行的矩函数{int i;for(i=1;i<=n+1;i++)A[0][i]=A[r][i];for(i=1;i<=n+1;i++)A[r][i]=A[k][i];for(i=1;i<=n+1;i++)A[k][i]=A[0][i];}float max(int k) //比校系数大小的函数{int i;float temp=0;for(i=k;i<=n;i++)if(fabs(A[i][k])>temp){temp=fabs(A[i][k]);flag=i;}return temp; }3.运行结果:截图如下:。
数值分析 迭代法 C++程序
课题三解线性方程组的迭代法实验目标:分别采用Jacobi迭代法,Gauss-Seidel迭代法和SOR迭代法求解线性方程组。
Jocabi迭代法:#include<iostream>#include<math.h>using namespace std;int i,j,k; //计数器int M = 2000;int Array(double ***Arr, int n){double **p;int i;p=(double **)malloc(n*sizeof(double *));if(!p)return 0;for(i=0;i<n;i++){p[i]=(double *)malloc(n*sizeof(double));if(!p[i])return 0;}*Arr=p;return 1;}void main(){double eps ;cout<<"默认最多迭代次数为2000次"<<endl<<"迭代精度为:";cin>>eps;double **matrix;int n;cout<<"矩阵大小为:";cin>>n;double *X;X= new double[n];double *Y;Y= new double[n];double *G;G= new double[n];for(i=0;i<n;i++){Y[i]=0;}if(!Array(&matrix,n))cout<<"内存分配失败!";elsecout<<"请输入矩阵:"<<endl;for( i=0;i<n;i++){for( j=0;j<n;j++){cin>>matrix[i][j];}}cout<<"请输入右端项:"<<endl;double *B;B = new double[n];for(i=0;i<n;i++){cin>>B[i];}for (i = 0 ;i< n;i++){if (fabs(matrix[i][i]) < eps){cout <<"打印失败"<<endl;return;}double T = matrix[i][i];for ( j = 0 ; j< n;j++){matrix[i][j] = -matrix[i][j]/T;}matrix[i][i] = 0;G[i] = B[i]/T;}int counter = 0;while (counter < M){for (i = 0;i < n; i++){double temp = 0;for (j = 0;j<n; j++){temp = temp + matrix[i][j]*Y[j];}X[i] = G[i] + temp;}double temp = 0;for (i = 0 ;i< n ; i++){temp = temp + fabs(X[i] - Y[i]);}if (temp <= eps)break;else{for( i = 0; i < n ;i++){Y[i] = X[i];}}counter++;}cout << "迭代次数为:"<<counter<<"次。
数值分析实验程序C语言
Y=Y+r*y[i];}
printf("%lf\n",Y);
}
2牛顿插值:
#include<stdio.h>
#define N 6
void main()
{
int i,k,n=N-1;
float X,Y,x[N],y[N],c[N],b[N];
printf("请输入xi:\n");
for(i=0;i<N;i++)
int k,N;
double x0,x1;
printf("请输入x0,N\n");
scanf("%lf,%d",&x0,&N);
for(k=1;k<=N;k++)
{
if(fd(x0)==0){printf("奇异标志\n");break;}
x1=x0-(f(x0)/fd(x0));
printf("X%d= %lf\n",k,x1);
}
}
float f(float x,float y)
{
float s;
s=y-(2*x)/y;
return s;
}
5牛顿迭代法:
#include<stdio.h>
#include<math.h>
#define e 1e-7
void main()
{
double f(double x);
double fd(double x);
1拉格朗日插值:
#include<stdio.h>
#define N 3
数值分析算法C语言程序
数值分析算法C语言程序数值分析是研究数学问题的近似解法的一门学科,其中包括了各种数值方法和算法。
本文将介绍数值分析中的常见算法,并给出相应的C语言程序。
1.二分法(Bisection Method)二分法是一种求函数零点的简单且常用的方法。
该方法的基本思想是通过不断将区间进行二分,并比较中点处函数值的正负来找到零点所在的区间。
```c#include <stdio.h>double f(double x)return x * x - 2;double bisection(double a, double b, double eps)double c;while ((b - a) > eps)c=(a+b)/2;if (f(c) == 0)break;}else if (f(a) * f(c) < 0)b=c;}elsea=c;}}return c;int maidouble a = 0.0;double b = 2.0;double result = bisection(a, b, eps);printf("The root is: %lf\n", result);return 0;```2.牛顿迭代法(Newton's Method)牛顿迭代法是一种高效的求函数零点的方法。
该方法的基本思想是通过对函数进行线性逼近,不断逼近函数的零点。
```c#include <stdio.h>#include <math.h>double f(double x)return x * x - 2;double df(double x)return 2 * x;double newton(double x0, double eps) double x = x0;double deltaX = f(x) / df(x);while (fabs(deltaX) > eps)deltaX = f(x) / df(x);x = x - deltaX;}return x;int maidouble x0 = 2.0;double result = newton(x0, eps); printf("The root is: %lf\n", result); return 0;```3.高斯消元法(Gaussian Elimination)高斯消元法是一种用于求解线性方程组的方法。
计算方法数值分析C语言源程序
第1章线性方程组的直接算法求解线性方程组的直接算法是基于矩阵分解的算法。
常见的矩阵分解有两种:1.矩阵的三角分解矩阵的三角就是把一个矩阵分解成两个三角形矩阵的乘积。
比如:简单的三角分解:,这里,是单位下三角矩阵,是上三角矩阵。
列主元三角分解:,这里,是初等置换阵,是单位下三角矩阵且各元素的模不超过1,是上三角矩阵。
全主元三角分解:,这里,,是初等置换阵,是单位下三角矩阵且元素的模不超过,是上三角矩阵。
2.矩阵的正交三角分解正交化三角化就是把一个矩阵分解成一个正交矩阵和一个上三角矩阵的乘积.即,这里是正交矩阵,是上三角矩阵.1.1 矩阵的三角分解1.1.1 功能把实矩阵分解成单位下三角形矩阵和上三角形矩阵的乘积.即.该算法适用于各阶顺序主子式不等于的矩阵.1.1.2 算法概述所谓三角分解就是把阶方阵作如下分解其中是单位下三角矩阵,上三角矩阵。
当时,构造Gauss变换则以此类推,只要对角线上的元素,便可以一直这样做下去。
直到将其化为上三角矩阵为止。
即有.其中,且从而有而且.综上所述,我们可以将两个矩阵因子继续存储在原矩阵的存储空间上.的主对角线上的1不予存储.算法5.3(计算三角分解:Gauss消去法)1.1.3 算法程序1.1.3.1.1 参数说明**a n阶矩阵; n 矩阵的阶;1.1.3.1.2 C程序bool GaussLU(double **a,int n)//n阶矩阵的LU分解{for(int k=0;k<n-1;k++){if(a[k][k]==0)return false;for(int i=k+1;i<n;i++){a[i][k]/=a[k][k];for(int j=k+1;j<n;j++)a[i][j]-=a[i][k]*a[k][j];}}return true;}1.1.4 例题求矩阵1 2 3 41 4 9 161 8 27 64的三角分解,结果如下:1 2 3 41 2 6 121 3 6 241 7 6 241.2 列主元三角分解1.2.1 功能用矩阵的列主元三角分解,分解矩阵:,这里,是初等置换阵,是单位下三角矩阵且各元素的模不超过1,是上三角矩阵。
数值分析中的各种公式C 代码
二分法2.2 算法步骤步骤1:准备计算f(x)在有根区间[a,b]端点处的值f(a),f(b).步骤2:二分计算f(x)在区间中点(a+b)/2处的值f((a+b)/2).步骤3:判断若f((a+b)/2)=0,则(a+b)/2即是根,计算过程结束,否则检验;若f((a+b)/2)f(a)<0,则以(a+b)/2代替b,否则以(a+b)/2代替a.反复执行步骤2和步骤3,直到区间[a,b]的长度小于允许误差e,此时中点(a+b)/2即为所求近似根。
2.3 程序流程图3 实验结果分析#include"stdio.h"void main(){float a,b,e,x,c;int k=0,n=1;scanf("%f,%f,%f",&a,&b,&e);while(n==1){x=(a+b)/2;c=(x*x*x-x-1)*(a*a*a-a-1);if(c<0){b=x;if(b-a<=e){printf("%f,%d\n",x,k);break;}else k++;}else{ if(c==0) { printf("%f,%d\n",x,k);break;}else { a=x; if(b-a<=e) { printf("%f,%d\n",x,k);break;}else k++;}}}}高斯塞德尔迭代法求方程组解 高斯主元素消去法求方程组解2.2 算法步骤高斯塞德尔迭代法:步骤1:给定初值)0(1)0(2)0(1,...,,n x x x ,精度ε,最大容许迭代次数M,令k=1。
步骤2:对 i=1,2,…,n 依此计算)0()1()0()1(01)0()1()1().(i i ii i iinj j j ij iix x x x e a x a xx→-=-=∑≠=步骤3:求出 e=}{max 1i ni e ≤≤,若 e<ε,则输出 )0(i x (i=1,2,..,n ),停止计算。
(整理)数值分析计算方法程序汇总
(一)秦九韶法例:已知p=2x2-3x-4,求x=1时,p=? 程序:#include "stdio.h"void main(){float a[50],b,x;int n,i,k;scanf("%d%f",&n,&x);for(i=0;i<=n;i++)scanf("%f",&a[i]);b=a[0];k=1;while(k<=n){b=x*b+a[k];k=k+1;}printf("%f",b);}结果:-5.000000(二)复化中矩形公式例:求() g y dy+∞⎰=[g(y1)+g(y2)+···+g(y m)]h的值(已知g(y)=y n/2-1e-y/2 ,n=11,h=0.1,y m=20)。
程序:#include "stdio.h"#include <math.h>main(){double y,h,s;int n,i;n=11;h=0.1;s=0;y=0;for(i=0;i<=200;i++){y=y+h;s=s+pow(y,n/2.0-1)*exp(-y/2.0);}s=s*h;printf("%f\n",s);}结果:2266.139467例:求121?x dx -=⎰程序:#include "stdio.h"void main(){double a,b,s,h,x;int i,n;a=-1.0;b=1.0;n=10;h=(b-a)/n;x=a;s=x*x/2;for(i=1;i<n;i++){x=x+h;s=s+x*x;}s=s+b*b/2;s=s*h;printf("s=%f\n",s);}结果:s=0.680000(四)复化辛普森公式例:求130?x dx=⎰程序:#include "stdio.h"void main(){double a,b,c,s,h,x,y;int i,n;a=0.0;b=1.0;n=10;s=0.0;h=(b-a)/n;x=a;y=x+h;c=(x+y)/2;for(i=1;i<=n;i++){s=s+x*x*x+4*c*c*c+y*y*y;x=x+h;y=y+h;c=c+h;}s=s*h/6;printf("s=%f\n",s);}结果:s=0.250000例:求220?x dx =⎰ 程序:#include <stdio.h>#include <math.h>main(){double a,b,h,s,x1,x2;int i,n;a=0;b=2;n=20;s=0;h=(b-a)/n;for(i=0;i<n;i++){x1=a+i*h+h/2*(1/1.732+1);x2=a+i*h+h/2*(1-1/1.732);s=s+x1*x1*x1+x2*x2*x2;}s=h/2*s;printf("s=%f\n",s);}结果:s=4.000000(六)二维中矩形公式例:求112200()?dx x y dy +=⎰⎰ 程序:#include <stdio.h>#include <math.h>main(){double a,b,c,d,x[50],y[50],hx,hy,s;int i,j,n,k;a=0.0;b=1.0;c=0.0;d=1.0;n=k=10;hx=(b-a)/n;hy=(d-c)/k;x[0]=a+hx/2;y[0]=c+hy/2;s=0.0;for(j=0;j<k;j++){for(i=0;i<n;i++){s=s+x[i]*x[i]+y[j]*y[j];x[i+1]=x[i]+hx;}y[j+1]=y[j]+hy;}s=hx*hy*s;printf("s=%f\n",s);}结果:s=0.665000(七)迭代法例:求x=x2的解。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
一、拉格朗日插值#include<stdio.h>#include<stdlib.h>#include<math.h>void Lagrange(float s){double x[5]={0.2,0.4,0.6,0.8,1.0},y[5]={0.98,0.92,0.81,0.64,0.38},f,L=0;int i,j;for (i=0;i<5;i++){ f=1; for (j=0;j<5;j++) if(j!=i) f=(s-x[j])/(x[i]-x[j])*f; L+=f*y[i]; } printf("输出:%f\n",L);}void main(){float x;printf("输入插值点:");scanf("%f",&x);Lagrange(x);}二、牛顿插值#include<stdlib.h>#include<stdio.h>#include<math.h>int ND(float s){double x[5]={0.2,0.4,0.6,0.8,1.0},y[5]={0.98,0.92,0.81,0.64,0.38},p=0,g,f;int i,j,k;for (i=0;i<5;i++){for (j=4;j>i;j--){ f=x[j]-x[j-i-1];y[j]=(y[j]-y[j-1])/f;}g=y[i+1];for (k=0;k<=i;k++)g=g*(s-x[k]);p=p+g;}printf("输出插值点函数值:%f\n",p+y[0]);return 1;}void main(){float x;printf("输入插值点:");scanf("%f",&x);ND(x);}三、埃尔米特插值#include<stdio.h>#include<stdlib.h>#include<math.h>void Hermite(float s){double x[3]={0.25,1,2.25},y[3]={0.125,1,3.375},z[3]={0.75,1.5,2.25};double H=0.0,a,b,f,g;int i,j;for (i=0;i<3;i++){f=1.0;g=0.0;for (j=0;j<3&&j!=i;j++){f=f*(s-x[j])/(x[i]-x[j]);g=g+1/(x[i]-x[j]);}a=(1-2*(s-x[i])*g)*f*f;b=(s-x[i])*f*f;H=H+y[i]*a+z[i]*b;}printf("%f\n",H);}void main(){float x;printf("输入插值点:");scanf("%f",&x);Hermite(x);}四、三次样条插值#include <math.h>#include <stdio.h>#include <stdlib.h>void main(){int N=7,R=2,i,k;double p1,p2,p3,p4;double x[8]={0.5,0.7,0.9,1.1,1.3,1.5,1.7,1.9};double y[8]={0.4794,0.6442,0.7833,0.8912,0.9636,0.9975,0.9917,0.9463};double P0=-0.4794,Pn=-0.9463,u[3]={0.6,0.8,1.2},s[3];double h[7],a[8],c[7], g[8],af[8],ba[7],m[8];for(k=0;k <N;k++)h[k]=x[k+1]-x[k];for(k=1;k <N;k++) a[k]=h[k]/(h[k]+h[k-1]);for(k=1;k <N;k++) c[k]=1-a[k];for(k=1;k <N;k++) g[k]=3*(c[k]*(y[k+1]-y[k])/h[k]+a[k]*(y[k]-y[k-1])/h[k-1]);c[0]=a[N]=1;g[0]=3*(y[1]-y[0])/h[0]-P0*h[0]/2;g[N]=3*(y[N]-y[N-1])/h[N-1]+Pn*h[N-1]/2;ba[0]=c[0]/2; g[0]=g[0]/2;for(i=1;i <N;i++){ af[i]=2-a[i]*ba[i-1]; g[i]=(g[i]-a[i]*g[i-1])/af[i]; ba[i]=c[i]/af[i]; }af[N]=2-a[N]*ba[N-1];g[N]=(g[N]-a[N]*g[N-1])/af[N];m[N]=g[N];for(i=N-1;i>=0;i--) m[i]=g[i]-ba[i]*m[i+1];for(i=0;i <=R;i++){ k=0;while(u[i]> x[k+1])k++;p1=(h[k]+2*(u[i]-x[k])*pow((u[i]-x[k+1]),2)*y[k])/pow(h[k],3);p2=(h[k]-2*(u[i]-x[k+1])*pow((u[i]-x[k]),2)*y[k+1])/pow(h[k],3);p3=(u[i]-x[k])*pow((u[i]-x[k+1]),2)*m[k]/pow(h[k],2);p4=(u[i]-x[k+1])*pow((u[i]-x[k]),2)*m[k+1]/pow(h[k],2);s[i]=p1+p2+p3+p4;}printf( "\nx= ");for(i=0;i <=N;i++)printf( "%8.1f ",x[i]);printf( "\ny= ");for(i=0;i <=N;i++)printf( "%8.4f ",y[i]);printf( "\n\nu= ");for(i=0;i <=R;i++)printf( "%9.2f ",u[i]);printf( "\n插值点:s= ");for(i=0;i <=R;i++)printf( "%9.5f ",s[i]);printf("\n");}五、复合梯形公式#include<stdio.h>#include<stdlib.h>#include<math.h>double FTX(int n,float a,float b){double f=0,t,h,*x,*y;int i;x=(double*)malloc((n+1)*sizeof(double));y=(double*)malloc((n+1)*sizeof(double));h=(b-a)/n;for(i=0;i<n+1;i++){x[i]=a+i*h;y[i]=sin(x[i]);}for(i=1;i<n;i++) f=f+2*y[i];t=h/2*(y[0]+f+y[n]);printf("输出函数值:%f\n",t);return 1;}void main(){float a,b;int n;printf("输入区间上,下限:");scanf("%f %f",&a,&b);printf("输入等分区间数:");scanf("%d",&n);FTX(n,a,b);}六、复合辛普森求积公式#include<stdio.h>#include<stdlib.h>#include<math.h>double FSP(int n,float a,float b){double f1=0,f2=0,h,*x1,*y1,*x2,*y2;int i;x1=(double*)malloc((n+1)*sizeof(double));y1=(double*)malloc((n+1)*sizeof(double));x2=(double*)malloc(n*sizeof(double));y2=(double*)malloc(n*sizeof(double));h=(b-a)/n;for(i=0;i<n+1;i++){x1[i]=a+h*i;y1[i]=sin(x1[i])*x1[i];}for(i=0;i<n;i++){x2[i]=x1[i]+h/2;y2[i]=sin(x2[i])*x2[i];}for(i=1;i<n;i++) f1=f1+2*y1[i];for(i=0;i<n;i++) f2=f2+4*y2[i];printf("输出函数值:%f\n",h/6*(y1[0]+f1+f2+y1[n]));return 1;}void main(){float a,b;int n;printf("输入区间上,下限:");scanf("%f %f",&a,&b);printf("输入等分区间数:");scanf("%d",&n);FSP(n,a,b);}七、直接三角分解法#include<stdio.h>#include<stdlib.h>#include<math.h>void main(){doubleA[3][3]={0.25,0.2,0.166667,0.3333,0.25,0.2,0.5,1,2},x[3],y[3],b[3]={9,8,8},L[3][3],U[3][3],f1=0,f2=0;int i,j,k;for(i=0;i<3;i++)for(j=0;j<3;j++){U[i][j]=0;L[i][j]=0;}for(i=0;i<3;i++){U[0][i]=A[0][i];L[i][0]=A[i][0]/U[0][0];L[i][i]=1;}for(i=1;i<3;i++)for(j=i;j<3;j++){for(k=0;k<=i-1;k++){f1=f1+L[i][k]*U[k][j];f2+=L[j][k]*U[k][i];} U[i][j]=A[i][j]-f1;L[j][i]=(A[j][i]-f2)/U[i][i];f1=0;f2=0;}y[0]=b[0];for(i=1;i<3;i++){for(j=0;j<=i-1;j++) f1+=L[i][j]*y[j];y[i]=b[i]-f1;f1=0;}x[2]=y[2]/U[2][2];for(i=1;i>=0;i--){for(j=i+1;j<3;j++) f2+=U[i][j]*x[j];x[i]=(y[i]-f2)/U[i][i];f2=0;} printf("输出L矩阵:\n");for(i=0;i<3;i++){for(j=0;j<3;j++)printf("%f ",L[i][j]);printf("\n");}printf("输出U矩阵:\n");for(i=0;i<3;i++){for(j=0;j<3;j++)printf("%f ",U[i][j]);printf("\n");}printf("输出求解结果:\n");for(i=0;i<3;i++)printf("%f ",x[i]);printf("\n");}八、改进的平方法#include<stdio.h>#include<stdlib.h>#include<math.h>void main(){double A[3][3]={2,-1,1,-1,-2,3,1,3,1},x[3],y[3],b[3]={4,5,6},d[3],L[3][3],U[3][3],f1=0,f2=0; int i,j,k,n=3;for(i=0;i<n;i++)for(j=0;j<n;j++) {U[i][j]=0;L[i][j]=0;}d[0]=A[0][0];L[0][0]=1;for(i=1;i<n;i++){L[i][i]=1;for(j=0;j<=i-1;j++) {for(k=0;k<=j-1;k++) f1+=U[i][k]*L[j][k];U[i][j]=A[i][j]-f1;L[i][j]=U[i][j]/d[j];f2+=U[i][j]*L[i][j];f1=0;}d[i]=A[i][j]-f2;f2=0;}y[0]=b[0];for(i=1;i<n;i++) {for(j=0;j<=i-1;j++) f1+=L[i][j]*y[j];y[i]=b[i]-f1;f1=0;}x[n-1]=y[n-1]/d[n-1];for(i=n-2;i>=0;i--) {for(j=i+1;j<n;j++) f2+=L[j][i]*x[j];x[i]=y[i]/d[i]-f2;f2=0;}printf("输出L矩阵:\n");for(i=0;i<n;i++){for(j=0;j<n;j++) printf("%f ",L[i][j]); printf("\n");}printf("输出U矩阵:\n");for(i=0;i<n;i++){for(j=0;j<n;j++) printf("%f ",U[i][j]); printf("\n");}printf("输出求解结果:\n");for(i=0;i<n;i++) printf("%f ",x[i]);printf("\n");}九、追赶法#include<stdio.h>#include<stdlib.h>#include<math.h>void main(){double A[5][5]={2,-1,0,0,0,-1,2,-1,0,0,0,-1,2,-1,0,0,0,-1,2,-1,0,0,0,-1,2},f[5]={1,0,0,0,0}; double x[5],y[5],U[5][5];int i,j,n=5;for(i=0;i<n;i++)for(j=0;j<n;j++) U[i][j]=0;U[0][0]=U[n-1][n-1]=1;U[0][1]=A[0][1]/A[0][0];for(i=1;i<n-1;i++){ U[i][i]=1;U[i][i+1]=A[i][i+1]/(A[i][i]-A[i][i-1]*U[i-1][i]);}y[0]=f[0]/A[0][0];for(i=1;i<n;i++)y[i]=(f[i]-A[i][i-1]*y[i-1])/(A[i][i]-A[i][i-1]*U[i-1][i]);x[n-1]=y[n-1];for(i=n-2;i>=0;i--) x[i]=y[i]-U[i][i+1]*x[i+1];printf("输出U矩阵:\n");for(i=0;i<n;i++){ for(j=0;j<n;j++) printf("%f ",U[i][j]); printf("\n");}printf("输出求解结果:\n");for(i=0;i<n;i++) printf("%f ",x[i]);printf("\n");}十、雅可比迭代法#include<stdio.h>#include<stdlib.h>#include<math.h>void main(){double A[3][3]={5,2,1,-1,4,2,2,-3,10},x[50][3],b[3]={-12,20,3},f=0,L=1;int n=3,i,j,k=1;for (i=0;i<3;i++) x[0][i]=0;/*初始值*/while(L<0.003){for(i=0;i<n;i++){ for(j=0;j<n;j++)if(j!=i) f=f+A[i][j]*x[k-1][j];x[k][i]=(b[i]-f)/A[i][i];}L=x[k][0]-x[k-1][0];for(i=1;i<n;i++)if((x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])>1||(x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])<-1) L=x[k][i]-x[k-1][i];k++;}printf("输出每次迭代求得的X值:\n");for(i=1;i<k;i++){printf("第%d次迭代:",i);for(j=0;j<n;j++)printf("%f ",x[i][j]);printf("\n");}printf("\n输出迭代次数:%d\n",k-1);}十一、高斯—塞德尔迭代法#include<stdio.h>#include<stdlib.h>#include<math.h>void main(){double A[3][3]={5,2,1,-1,4,2,2,-3,10},x[50][3],b[3]={-12,20,3},f1=0,f2=0,L=1;int i,j,k=1,n=3;for (i=0;i<3;i++) x[0][i]=0;while(L>0.00001||L<-0.00001){for(i=0;i<n;i++){ for(j=0;j<n;j++){if(j<i) f1+=A[i][j]*x[k][j];else if(j>i) f2+=A[i][j]*x[k-1][j];}x[k][i]=(b[i]-f1-f2)/A[i][i];}L=x[k][0]-x[k-1][0];for(i=1;i<n;i++)if((x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])>1||(x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])<-1) L=x[k][i]-x[k-1][i];k++;}printf("输出迭代X矩阵:\n");for(i=1;i<k;i++){for(j=0;j<n;j++) printf("%f ",x[i][j]); printf("\n");}printf("\n输出迭代次数:%d\n",k-1);}十二、超松弛迭代法#include<stdio.h>#include<stdlib.h>#include<math.h>void main(){double A[3][3]={5,2,1,-1,4,2,2,-3,10},x[30][3],b[3]={-12,20,3},f1=0,f2=0,L=1,w=0.9;int i,j,k=1,n=3;for(i=0;i<3;i++) x[0][i]=0;while(L>0.0001||L<-0.0001){for(i=0;i<n;i++){ for(j=0;j<n;j++){if(j<i) f1+=A[i][j]*x[k][j];else if(j>i) f2+=A[i][j]*x[k-1][j];}x[k][i]=w*(b[i]-f1-f2)/A[i][i];}L=x[k][0]-x[k-1][0];for(i=1;i<n;i++)if((x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])>1||(x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])<-1) L=x[k][i]-x[k-1][i];k++;}printf("输出迭代X矩阵:\n");for(i=1;i<k;i++){for(j=0;j<n;j++) printf("%f ",x[i][j]); printf("\n");}printf("\n输出迭代次数:%d\n",k-1);}数值分析算法程序班级:计算08Q2班姓名:甄彦福。