/******************************************************************************
Online C Compiler.
Code, Compile, Run and Debug C program online.
Write your code in this editor and press "Run" button to compile and execute it.
*******************************************************************************/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#define PI 3.141592
void read_input(double *D, double *L, int *nx, double *t_F);
double main(void) {
/******************************/
/* Declarations of parameters */
/******************************/
/* Number of grid points */
int nx;
/* Length of domain */
double L;
/* Equation coefficients */
double D;
/* Length of time to run simulation. */
double t_F;
/* Read in from file; */
read_input(&D, &L, &nx, &t_F);
/* Grid spacing */
double dx = L/nx;
double invdx2 = 1.0/(dx*dx);
/* Time step */
double dt = 0.25/invdx2; // changed to 0.25/dx^2 to satisfy the stability condition
/************************************************/
/* Solution Storage at Current / Next time step */
/************************************************/
double *uc, *un, *vc, *vn;
/* Time splitting solutions */
double *uts1, *uts2, *vts1, *vts2;
/* Derivative used in finite difference */
double deriv;
/* Allocate memory according to size of nx */
uc = malloc(nx * sizeof(double));
un = malloc(nx * sizeof(double));
vc = malloc(nx * sizeof(double));
vn = malloc(nx * sizeof(double));
uts1 = malloc(nx * sizeof(double));
uts2 = malloc(nx * sizeof(double));
vts1 = malloc(nx * sizeof(double));
vts2 = malloc(nx * sizeof(double));
/* Check the allocation pointers */
if (uc==NULL||un==NULL||vc==NULL||vn==NULL||uts1==NULL||
uts2==NULL||vts1==NULL||vts2==NULL) {
printf("Memory allocation failed\n");
return 1;
}
int k;
double x;
/* Current time */
double ctime;
/* Initialise arrays */
for(k = 0; k < nx; k++) {
x = k*dx;
uc[k] = 1.0 + sin(2.0*PI*x/L);
vc[k] = 0.0;
/* Set other arrays to 0 */
uts1[k] = 0; uts2[k] = 0;
vts1[k] = 0; vts2[k] = 0;
}
/* Loop over timesteps */
while (ctime < t_F){
/* Rotation factors for time-splitting scheme. */
double cfac = cos(dt); //changed from 2*dt to dt
double sfac = sin(dt);
/* First substep for diffusion equation, A_1 */
for (k = 0; k < nx; k++) {
x = k*dx;
/* Diffusion at half time step. */
deriv = (uc[k-1] + uc[k+1] - 2*uc[k])*invdx2 ;
uts1[k] = uc[k] + (D * deriv + vc[k])* 0.5*dt; //
deriv = (vc[k-1] + vc[k+1] - 2*vc[k])*invdx2;
vts1[k] = vc[k] + (D * deriv - uc[k]) * 0.5*dt;
}
/* Second substep for decay/growth terms, A_2 */
for (k = 0; k < nx; k++) {
x = k*dx;
/* Apply rotation matrix to u and v, */
uts2[k] = cfac*uts1[k] + sfac*vts1[k];
vts2[k] = -sfac*uts1[k] + cfac*vts1[k];
}
/* Third substep for diffusion terms, A_1 */
for (k = 0; k < nx; k++) {
x = k*dx;
deriv = (uts2[k-1] + uts2[k+1] - 2*uts2[k])*invdx2;
un[k] = uts2[k] + (D * deriv + vts2[k]) * 0.5*dt;
deriv = (vts2[k-1] + vts2[k+1] - 2*vts2[k])*invdx2;
vn[k] = vts2[k] + (D * deriv - uts2[k]) * 0.5*dt;
}
/* Copy next values at timestep to u, v arrays. */
memcpy(uc,un, sizeof(double) * nx);
memcpy(vc,vn, sizeof(double) * nx);
/* Increment time. */
ctime += dt;
for (k = 0; k < nx; k++ ) {
x = k*dx;
printf("%g %g %g %g\n",ctime,x,uc[k],vc[k]);
}
}
/* Free allocated memory */
free(uc); free(un);
free(vc); free(vn);
free(uts1); free(uts2);
free(vts1); free(vts2);
return 0;
}
// The lines below don't contain any bugs! Don't modify them
void read_input(double *D, double *L, int *nx, double *t_F) {
FILE *infile;
if(!(infile=fopen("input.txt","r"))) {
printf("Error opening file\n");
exit(1);
}
if(4!=fscanf(infile,"%lf %lf %d %lf",D,L,nx,t_F)) {
printf("Error reading parameters from file\n");
exit(1);
}
fclose(infile);
}
所以这是代码。它旨在解决以下微分方程：
du/dt - Dd^2u/dx^2 - v = 0
dv/dt - Dd^2v/dx^2 + u = 0
它把方程分成两部分。第二个 x 导数部分(A1) 和包含 u 和 v(A2) 的衰减部分。它对 A1 使用两个半步（0.5dt），对 A2 使用 1 个全步 dt。我知道如何进行时间分割，但我不知道我是否在这里正确地完成了。
这是一个作业，我已经修复了所有错误，我只是想让代码按预期工作。我从来没有解决过类似的问题，所以我现在肯定非常卡住。解决方案收敛，但我认为它是错误的。任何想法为什么？我不是在寻找直接告诉我做错了什么的人，如果您知道我的意思，请引导我朝着正确的方向前进。
PS：当我用 gcc 编译代码时，我收到一个关于双主（void）的警告。为什么会这样？