#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;}。
本科生课程设计报告实习课程数值分析学院名称管理科学学院专业名称信息与计算科学学生姓名学生学号指导教师实验地点实验成绩二〇一六年六月二〇一六年六摘要常微分方程数值解法是计算数学的一个分支.是解常微分方程各类定解问题的数值方法.现有的解析方法只能用于求解一些特殊类型的定解问题，实用上许多很有价值的常微分方程的解不能用初等函数来表示，常常需要求其数值解.所谓数值解，是指在求解区间内一系列离散点处给出真解的近似值.这就促成了数值方法的产生与发展.?关键词：数值解法；常微分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 预测-校正系统，利用外插公式计算预测，用内插公式进行校正。
#include <iostream>#include<stdio.h>#include<stdlib.h>#include<malloc.h>#include<string.h>using namespace std;#define OK 1#define ERROR 0#define OVERFLOW -2#define STACK_SIZE 100#define STACKINCREMENT 10#define OPSIZE 7typedef double ElemType;typedef char SElemType;char op[OPSIZE]={'+','-','*','/','(',')','#'};charprior[OPSIZE][OPSIZE]={{'>','>','<','<','<','>','>'},{'>','>','<','<','<','>','>'},{'>','>','>','>','<','>','>'},{'>','>',' >','>','<','>','>'},{'<','<','<','<','<','=',' '},{'>','>','>','>',' ','>','>'},{'<','<','<','<','<',' ','='}};typedef struct {ElemType *top;ElemType *base;int Stacksize;}OPND;typedef struct {SElemType *top;SElemType *base;int Stacksize;}OPTR;void InitStack(OPND *S1){S1->base=(ElemType *)malloc(STACK_SIZE * sizeof(ElemType));if(NULL==S1->base){printf("Init Failed!\n");exit (OVERFLOW);}else{S1->top=S1->base;S1->Stacksize=STACK_SIZE;}}int StackEmpty(OPND *S1){if(S1->top==S1->base)return 1;elsereturn 0;}int StackEmpty(OPTR *S2){if(S2->top==S2->base)return 1;elsereturn 0;}void InitStack(OPTR *S2){S2->base=(SElemType *)malloc(STACK_SIZE*sizeof(SElemType));if(NULL==S2->base){printf("Init Failed!\n");exit (OVERFLOW);}else{S2->top=S2->base;S2->Stacksize=STACK_SIZE;}}void Push(OPND *S1,ElemType e){if(S1->top-S1->base>=S1->Stacksize){S1->base=(ElemType*)realloc(S1->base,(S1->Stacksize+STACKINCREMENT)*sizeof(ElemType));if(!S1->base)exit(OVERFLOW);S1->top=S1->base+S1->Stacksize;S1->Stacksize=S1->Stacksize+STACKINCREMENT;}*S1->top=e;S1->top++;}void Push(OPTR *S2,SElemType e){if(S2->top-S2->base>=S2->Stacksize){S2->base=(SElemType*)realloc(S2->base,(S2->Stacksize+STACKINCREMENT)*sizeof(SElemType));if(!S2->base)exit (OVERFLOW);S2->top=S2->base+S2->Stacksize;S2->Stacksize+=STACKINCREMENT;}*S2->top++=e;}ElemType Pop(OPND *S1,ElemType *e){if(!StackEmpty(S1)){S1->top=S1->top-1;*e=*S1->top;}else{return ERROR;}return OK;}SElemType Pop(OPTR *S2,SElemType *e){if(StackEmpty(S2)){return ERROR;}else{*e=*--S2->top;}return OK;}ElemType GetTop(OPND *S1,ElemType *e){if(!StackEmpty(S1)){*e=*(S1->top-1);}return OK;}SElemType GetTop(OPTR *S2,SElemType *e){if(!StackEmpty(S2)){*e=*(S2->top-1);}return OK;}bool Judege(char c,char *op)//判断C是否为运算符{bool check=false;int i;for(i=0;i<OPSIZE;i++){if(c==op[i]){check = true;}}return check;}int ReturnOpOrd(char c,char *op)//返回运算符c在OP中的位置{int i;for(i=0;i<OPSIZE;i++){if(c==op[i]){return i;}}return OK;}char Precede(char theat1,char theat2)//g判断运算符栈的栈顶运算符theta1和读入的运算符theta2优先关系{return prior[ReturnOpOrd(theat1,op)][ReturnOpOrd(theat2,op)];}double Operate(ElemType a,SElemType theat,ElemType b)//操作数进行运算{switch(theat){case '+':return a+b;break;case '-':return a-b;break;case '*':return a*b;break;case '/':if(b!=0){return a/b;}else{printf("出错！除数为0！\n");return ERROR;}break;default:return 0;}}float EvaluateExpression(char *expression){OPND S1;OPTR S2;char *c=expression;ElemType e1,e2,result,d;char temp[10],r[2];SElemType e;InitStack(&S1);InitStack(&S2);Push(&S2,'#');strcpy(temp,"\0");while(*c!='#'||((GetTop(&S2,&e),e))!='#'){if(!Judege(*c,op)){r[0]=*c;r[1]='\0';strcat(temp,r);c++;if(Judege(*c,op)){d=(float)atof(temp);Push(&S1,d);strcpy(temp,"\0");}}else{switch(Precede( (GetTop(&S2,&e),e),*c)){case '<':Push(&S2,*c);c++;break;case '=':Pop(&S2,&e);c++;break;case '>':Pop(&S1,&e1);Pop(&S1,&e2);Pop(&S2,&e);result=Operate(e2,e,e1);Push(&S1,result);break;}}}Pop(&S1,&result);return result;}int main(){char expression[80];float result;printf("输入你要计算的算式：\n");gets(expression);result=EvaluateExpression(expression);printf("计算结果为：%f",result);return OK;}。
#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;2;t[0][0]=h*(f(a)+f(b))/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);2;t[0][k]=(t[0][k]+t[0][k-1])/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)得到近似结果为 x ≈ %lf\n\n",x,i);printf("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;初值为 (x0,y0) = ( %.8f , %.8f )\n",x,y);printf("\nfor(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);6;y=y+h*(k1+2*(k2+k3)+k4)/x=x+h;第%d次输出结果为 (x%d,y%d) = ( %.8f , %.8f )\n",i+1,i+1,i+1,x,y);printf("}}#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];}在 x = %f 处的近似值为: y = %f\n",x,y);printf("\n f(x)}#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];}方程组的解为:\n");printf("for(j=0;j<m+1;j++)printf("a%d = %f\n",j,a[j][m+1]);拟合多项式为:\n");printf("printf("P%d(x) = (%f) + (%f)x + (%f)x^2\n",m,a[0][m+1],a[1][m+1],a[2][m+1]);}else数据有误！\n");printf("}}列主元素法#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;}方程第%d行互换位置后如下\n",i+1);printf("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;}方程经%d次消元如下\n",i+1);printf("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){方程化简得\n");printf("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];}方程组的解为:\n");printf("for(j=0;j<3;j++)printf("x%d = %f\n",j+1,a[j][3]);}elseprintf("数据有误！\n");}#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;//把a矩阵转化为b矩阵//for(i=0;i<3;i++){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;}化为b矩阵如下\n");printf("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){计算结果为\n");printf("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);第%d次输出结果为(%.8f,%.8f)\n",n,x,y);printf("x0=x;y0=y;}}。
