You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
150 lines
3.3 KiB
150 lines
3.3 KiB
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <math.h>
|
|
typedef struct {
|
|
int N;
|
|
double d;
|
|
double h;
|
|
double xmin, xmax, tmin, tmax;
|
|
double j;
|
|
double n;
|
|
int jmax;
|
|
int nmax;
|
|
|
|
} soli_param;
|
|
double *u0;
|
|
double *uj1;
|
|
void ustart(soli_param *sp) {
|
|
double N = sp->N;
|
|
double jmax = sp->jmax;
|
|
double h = sp->h;
|
|
double x = sp->xmin;
|
|
|
|
for (int j = 0; j <= jmax; j++) {
|
|
|
|
u0[j] = -N * (N + 1) / (pow(cosh(x), 2));
|
|
x = x + h;
|
|
printf("%f\n", u0[j]);
|
|
|
|
}
|
|
|
|
}
|
|
void ujeins(soli_param *sp) {
|
|
double jmax = sp->jmax;
|
|
double d = sp->n;
|
|
double h = sp->h;
|
|
for (int j = 0; j <= jmax; j++) {
|
|
if (j == 1) {
|
|
uj1[j] = (d * 6 * u0[j] * (u0[j + 1] - u0[j - 1]) / (2 * h))
|
|
- (d * (u0[j + 2] + 2 * u0[j - 1]) / (2 * pow(h, 3)))
|
|
+ u0[j];
|
|
|
|
} else if (j == 0) {
|
|
uj1[j] = (d * 6 * u0[j] * (u0[j + 1] - u0[j - 1]) / (2 * h))
|
|
- (d * (u0[j + 2]) / (2 * pow(h, 3))) + u0[j];
|
|
|
|
} else if (j == jmax - 1) {
|
|
uj1[j] = (d * 6 * u0[j] * (u0[j + 1] - u0[j - 1]) / (2 * h))
|
|
- (d * (2 * u0[j - 1] - u0[j - 2]) / (2 * pow(h, 3)))
|
|
+ u0[j];
|
|
|
|
} else if (j == jmax) {
|
|
uj1[j] = (d * 6 * u0[j] * (-u0[j - 1]) / (2 * h))
|
|
- (d * (2 * u0[j - 1] - u0[j - 2]) / (2 * pow(h, 3)))
|
|
+ u0[j];
|
|
} else {
|
|
uj1[j] = (d * 6 * u0[j] * (u0[j + 1] - u0[j - 1]) / (2 * h))
|
|
- (d * (u0[j + 2] + 2 * u0[j - 1] - u0[j - 2])
|
|
/ (2 * pow(h, 3))) + u0[j];
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
void nextu(soli_param *sp) {
|
|
int jmax = sp->jmax;
|
|
int nmax = sp->nmax;
|
|
double d = sp->d;
|
|
double h = sp->h;
|
|
double * unm = u0;
|
|
double * un = uj1;
|
|
double unp[jmax];
|
|
for (int n = 0; n <= nmax; n++) {
|
|
for (int j = 0; j <= jmax; j++) {
|
|
if (j == 1) {
|
|
unp[j] = 2*d*(
|
|
(un[j + 1] + un[j] + un[j - 1])
|
|
* (un[j + 1] - un[j - 1]) / h
|
|
- (un[j + 2] - 2 * un[j + 1]
|
|
+ 2 * un[j - 1])
|
|
/ (2 * pow(h, 3))) + unm[j];
|
|
|
|
} else if (j == 0) {
|
|
unp[j] = 2
|
|
* d*(
|
|
(un[j + 1] + un[j]) * (un[j + 1]) / h
|
|
- (un[j + 2] - 2 * un[j + 1])
|
|
/ (2 * pow(h, 3))) + unm[j];
|
|
|
|
} else if (j == jmax - 1) {
|
|
unp[j] = 2
|
|
* d*(
|
|
(un[j + 1] + un[j] + un[j - 1])
|
|
* (un[j + 1] - un[j - 1]) / h
|
|
- (-2 * un[j + 1] + 2 * un[j - 1]
|
|
- un[j - 2]) / (2 * pow(h, 3)))
|
|
+ unm[j];
|
|
|
|
} else if (j == jmax) {
|
|
unp[j] = 2
|
|
* d*(
|
|
(un[j] + un[j - 1]) * (-un[j - 1]) / h
|
|
- (2 * un[j - 1] - un[j - 2])
|
|
/ (2 * pow(h, 3))) + unm[j];
|
|
|
|
} else {
|
|
unp[j] = 2
|
|
* d*(
|
|
(un[j + 1] + un[j] + un[j - 1])
|
|
* (un[j + 1] - un[j - 1]) / h
|
|
- (un[j + 2] - 2 * un[j + 1]
|
|
+ 2 * un[j - 1] - un[j - 2])
|
|
/ (2 * pow(h, 3))) + unm[j];
|
|
}
|
|
|
|
|
|
}
|
|
unm=un;
|
|
un=unp;
|
|
|
|
}
|
|
FILE *fp= fopen("numplot2.csv", "w");
|
|
double x=sp->xmin;
|
|
for(int j=0;j<=jmax;j++){
|
|
x=x+h;
|
|
fprintf(fp, "%f\t%f\n",x,un[j]);
|
|
printf( "%f\t%f\n",x,un[j]);
|
|
|
|
}
|
|
fclose(fp);
|
|
}
|
|
|
|
|
|
int main(void) {
|
|
soli_param sp;
|
|
sp.N = 2;
|
|
sp.h=0.1;
|
|
sp.d=pow(sp.h,3)/3;
|
|
sp.tmax=1;
|
|
sp.tmin=0;
|
|
sp.xmin=-30;
|
|
sp.xmax=30;
|
|
sp.jmax=(sp.xmax-sp.xmin)/sp.h;
|
|
sp.nmax=(sp.tmax-sp.tmin)/sp.d;
|
|
u0 = (double*) malloc(sp.jmax * sizeof(double));
|
|
ustart(&sp);
|
|
uj1 = (double*) malloc(sp.jmax * sizeof(double));
|
|
ujeins(&sp);
|
|
nextu(&sp);
|
|
|
|
}
|