gpt4 book ai didi

c++ - 在 C++ 中实现 crank-nicolson

转载 作者:行者123 更新时间:2023-11-27 23:28:53 24 4
gpt4 key购买 nike

我想根据那个在 C++ 中实现 crank-nicolson 方法:

−ru(i−1,j+1) + 2(1 + r)u(i,j+1) − ru(i+1,j+1) = 2(1 − r)u(i,j) + r (u(i−1,j)+ u(i+1,j) )

我使用 gauss-jordan 方法求解系统,但我不知道如何实现上述公式。

const double pi=3.14159265;

double f (double x){
return sin(pi*x);

}

using namespace std;


//gauss-jordan
double* gauss(int n ,double **a){
double factor;
double *b,*x;
x=new double[n];
b=new double[n];


for (int k=1;k<=n;k++) {
for (int i=1;i<=n;i++) {
if (i!=k) {
factor=a[i][k]/a[k][k];
for (int j=1; j<=n; j++ ) {
a[i][j]=a[k][j]*factor-a[i][j];

}
b[i]=b[k]*factor -b[i];
}

}
}
for (int i=1;i<=n;i++) {
x[i]=b[i]/a[i][i];

}
return x;
}

int main()
{
int i,j,n,m,xd,td;
double h,k,r;
double **t,**p;

//----- initialize double pointer------------------
double **u;
u=new double *[m];
for (int o=1;o<=m;o++){
u[o]=new double [n];
}
//----- end of initialization-----------------------

cout <<"\nEnter the value of x dimension : \n";
cin >> xd;
cout <<"\nEnter the step size h for x dimension : \n";
cin >>h;
cout <<"\nEnter the value of time dimension : \n";
cin >>td;
cout<<"\nEnter the step k of time dimension : \n";
cin >>k;
n=xd/h -1.0;
m=td/k -1.0;
cout <<"\n\nThe internal elements of x dimension are :"<<n;
cout <<"\nThe internal elements of t dimension are :"<<m;
r=k/(h*h);
cout <<"\nThe r value is : "<<r;


//initial conditions
for (j=0;j<=m;j++){
u[0][m]=0;
u[10][m]=0;
}

//get the function
for (i=1;i<n;i++){
u[i][0]=f(i*h);
}


//apply crank-nicolson
for (i=1;i<n;i++){
for (j=1;j<n;j++){
-r*u[i-1][j+1] +2.0*(1.0+r)*u[i][j+1] -r*u[i+1][j+1]=2.0*(1.0-r)*u[i][j] +r*(u[i-1][j]+u[i+1][j]);
} // here i can't figure the steps i must follow in order for this to work



//-----delete double pointer-------------
for(int o=1;o<m;o++){
delete [] u[o];
delete [] u;
}
//---------------------------------------


return 0;
}

最佳答案

我假设变量 j 代表时间步长。为了实现 Crank-Nicolson,您必须将问题作为线性方程组提出并求解。系统对应的矩阵会是三对角形式的,所以最好用Thomas' algorithm而不是 Gauss-Jordan。

线性系统的形式为 A x = b,其中 x 是 vector (..., u(i-1, j+1), u(i, j+1), u(i+ 1, j+1), ...) 和 b 是 vector (..., r u(i−1,j), 2(1 − r)u(i,j), r u(i+1,j ), ...)。矩阵 A 的第 i 行的形式为(0, ..., 0, -r, 2(1 + r), -r, 0, ..., 0)。您必须小心处理第一行和最后一行,您必须在此处替换问题的边界条件。

关于有限差分法,尤其是 Crank-Nicolson 的一个很好的引用是 book by John Strikwerda .

希望这对您有所帮助。

关于c++ - 在 C++ 中实现 crank-nicolson,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/7240380/

24 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com