微分方程数值解法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.youerw.com 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.youerw.com 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.youerw.com 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
