毕业论文论文范文课程设计实践报告法律论文英语论文教学论文医学论文农学论文艺术论文行政论文管理论文计算机安全
您现在的位置: 毕业论文 >> 课程设计 >> 正文

微分方程数值解法C语言

更新时间:2010-7-14:  来源:毕业论文
微分方程数值解法C语言
由于对matlab语言不熟悉,所以还是采用C。前面几个都比较简单,最后一个需要解非其次方程组。采用高斯—Jordan消元法(数值分析)求逆解方程组,也再一次体会到算法本身的重要性,而不是语言。当然,矩阵求逆的算法也在100个经典的C语言算法之列。不过偏微分方程数值解的内容的确比较高深,我只能停留在编这种低级的东西的自娱自乐中。不过解决计算机、数学、信计专业的课程设计还是足够了。由于篇幅所限,只把源代码粘贴在这。
  一:预报矫正格式
  #include <math.h>
  #include<iostream>
  #include<stdlib.h>
  double count_0( double xn,double yn){
  //矫正格式
  double s;
  s=yn+0.1*(yn/xn*0.5+xn*xn/yn*0.5);
  return s;
  }
  double count_1(double xn,double yn,double y0){
  //预报格式
  double s;
  s=yn+0.05*((yn/xn*0.5+xn*xn/yn*0.5)+(y0/xn*0.5+xn*xn/y0*0.5));
  return s;
  }
  void main(){
  //计算,步长为0.1,进行10次计算,设初始值毕业论文http://www.751com.cn
  double xn=1,yn=1;
  int i=1;
  while(i<=10){
  printf("%16f ,%1.16f ,%1.16f \n",xn,yn,count_1(xn,yn,count_0(xn,yn)));
  xn=xn+0.1;
  yn=count_1(xn,yn,count_0(xn,yn));
  i++;
  }
  }
  二 显示差分格式
  #include<iostream>
  #include<math.h>
  #include<stdlib.h>
  main(){
  double a[6][11];
  //初始化;
  for(int i=0;i<=5;i++)
  {a[0]=0;a[10]=0;}
  for(int j=1;j<10;j++)
  {
  double p=3.14*j*0.1;
  a[0][j]=sin(p);
  }
  //按显示格式计算
  for(i=1;i<=5;i++)
  for(j=1;j<10;j++)
  a[j]=a[i-1][j-1]+a[i-1][j+1];
  //输出计算好的矩阵
  for(i=0;i<=5;i++)
  {
  for(j=0;j<11;j++)
  printf("%1.10f ",a[j]);
  printf("\n");
  }
  }
  三 龙阁库塔格式
  #include <math.h>
  #include<iostream>
  #include<stdlib.h>
  double count_k( double xn,double yn){
  double s;
  s=yn/xn*0.5+xn*xn/yn*0.5;
  return s;
  }
  void main(){
  //步长为0.1
  double xn=1,yn=1;
  int i=1;
  while(i<=11){
  printf("%f ,%f\n",xn,yn);
  double k1=count_k(xn,yn);
  double k2=count_k(xn+0.05,yn+0.05*k1);
  double k3=count_k(xn+0.05,yn+0.05*k2);
  double k4=count_k(xn+0.01,yn+0.1*k3);
  yn=yn+0.1/6*(k1+2*k2+2*k3+k4);
  xn=xn+0.1;
  i++;
  }
  }
    
  四 CRANK--NICOLSON隐式格式
  #include<iostream>
  #include<math.h>
  #include<stdlib.h>
  double Surplus(double A[],int m,int n);
  double * MatrixInver(double A[],int m,int n);
  double * MatrixOpp(double A[],int m,int n) /*矩阵求逆*/
  {
  int i,j,x,y,k;
  double *SP=NULL,*AB=NULL,*B=NULL,X,*C;
  SP=(double *)malloc(m*n*sizeof(double));
  AB=(double *)malloc(m*n*sizeof(double));
  B=(double *)malloc(m*n*sizeof(double));
  X=Surplus(A,m,n);
  X=1/X; 毕业论文http://www.751com.cn
  for(i=0;i<m;i++)
  for(j=0;j<n;j++)
  {for(k=0;k<m*n;k++)
  B[k]=A[k];
  {for(x=0;x<n;x++)
  B[i*n+x]=0;
  for(y=0;y<m;y++)
  B[m*y+j]=0;
  B[i*n+j]=1;
  SP[i*n+j]=Surplus(B,m,n);
  AB[i*n+j]=X*SP[i*n+j];
  }
  }
  C=MatrixInver(AB,m,n);
  return C;
  }
  double * MatrixInver(double A[],int m,int n) /*矩阵转置*/
  {
  int i,j;
  double *B=NULL;
  B=(double *)malloc(m*n*sizeof(double));
  for(i=0;i<n;i++)
  for(j=0;j<m;j++)
  B[i*m+j]=A[j*n+i];
  return B;
  }
  double Surplus(double A[],int m,int n) /*求矩阵行列式*/
  {
  int i,j,k,p,r;
  double X,temp=1,temp1=1,s=0,s1=0;
  if(n==2)
  {for(i=0;i<m;i++)
  for(j=0;j<n;j++)
  if((i+j)%2) temp1*=A[i*n+j];
  else temp*=A[i*n+j];
  X=temp-temp1;}
  else{
  for(k=0;k<n;k++)
  {for(i=0,j=k;i<m,j<n;i++,j++)
  temp*=A[i*n+j];
  if(m-i)
  {for(p=m-i,r=m-1;p>0;p--,r--)
  temp*=A[r*n+p-1];}
  s+=temp;
  temp=1;
  }
  for(k=n-1;k>=0;k--)
  {for(i=0,j=k;i<m,j>=0;i++,j--)
  temp1*=A[i*n+j];
  if(m-i)
  {for(p=m-1,r=i;r<m;p--,r++)
  temp1*=A[r*n+p];}
  s1+=temp1;
  temp1=1;
  }
  X=s-s1;}
  return X;
  }
  void initmat_A(double a[][9],double r){
  for(int i=0;i<9;i++)
  for(int j=0;j<9;j++)
  a[j]=0;
  for(i=0;i<9;i++){
  a=1+r;
  if(i!=8) a[i+1]=-0.5*r;
  if(i!=0) a[i-1]=-0.5*r;
  }
  }
  void initmat_B(double b[][9],double r){
  for(int i=0;i<9;i++)
  for(int j=0;j<9;j++)
  b[j]=0;
  for( i=0;i<9;i++){
  b=1-r;
  if(i!=8) b[i+1]=0.5*r;
  if(i!=0) b[i-1]=0.5*r;
  }
  }
  void initmat_C(double C[][9]){
  for(int i=0;i<9;i++)
  for(int j=0;j<9;j++)
  C[j]=0;
  }
  void main(){
  double a[100][11];
  for(int i=0;i<100;i++)
  for(int j=0;j<11;j++)
  a[j]=0;
  //初始化;
  for(i=0;i<100;i++)
  {a[0]=0;a[10]=0;}
  for(int j=1;j<10;j++)
  {
  double p=4*3.14*j*0.1;
  a[0][j]=sin(p);
  }
  //取h=0.1*3.14,r=0.0005,t=0.0001*3.14*3.14;
  //得到矩阵a和矩阵b
  double A[9][9];initmat_A(A,0.005);
  double B[9][9];initmat_B(B,0.005);
  //B矩阵与Un相乘,en是0;
  double C[9][9];initmat_C(C);
  double *A_;
  A_=MatrixOpp(A[0],9,9);//A矩阵求逆;
  //A逆*B 毕业论文http://www.751com.cn
  for(i=0;i<9;i++)
  for(j=0;j<9;j++)
  for(int s=0;s<9;s++)
  C[j]+=A_[i*9+s]*B[s][j];
  //填写a表格
  for(i=0;i<100;i++)
  {
  for(j=1;j<10;j++)
  for(int s=0;s<9;s++)
  a[i+1][j]+=a[s+1]*C[j-1][s];
  }
  //输出表格
  for(i=0;i<100;i++)
  {
  for(j=0;j<11;j++)
  printf("%1.8f ",a[j]);
  printf("\n");
  }
  printf("\n"); printf("\n");
  //利用精确解,求出表格
  for(i=0;i<100;i++)
  {
  for(j=0;j<11;j++)
  printf("%1.8f ",exp(-16*0.0001*0.0005*3.14*3.14*i)*sin(4*j*0.1*3.14));
  printf("\n");
  }
  }
微分方程数值解法C语言下载如图片无法显示或论文不完整,请联系qq752018766
设为首页 | 联系站长 | 友情链接 | 网站地图 |

copyright©751com.cn 辣文论文网 严禁转载
如果本毕业论文网损害了您的利益或者侵犯了您的权利,请及时联系,我们一定会及时改正。