微分方程数值解法C语言-课程设计
微分方程数值解法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);printf输出格式matlab
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次计算，设初始值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;
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;