omm

#include <stdio.h>

#include <math.h>

#define n 500

#define pi 3.14159265

void solve(double);

double g(double);

double dg(double);

double f(double,double,double,double);

double df(double);

double ux0(double);

double ut0(double);

FILE *fp;

double xn=0.0,xk=-1.0; // нач. и кон. знач

double dt,h,eps=0.0001; // шаг по времени, по координате, точность,

double T=0.3, Nt=300; // конечное время, кол-во шагов по времени

double u[n+1], un[n+1]; // массивы значений u(x)

int i, k;

void main(void)

{

fp=fopen("data.txt","w");

dt=T/Nt; // шаг по времени

h=(xk-xn)/n; // шаг по х

for(i=0;i<=n;i++)

{

u[i]=ut0(h*i); // нач. усл.

if(i%10==0) fprintf(fp,"%f\t%f\t%f\n",0.0,i*h,u[i]); // печать в файл начальных значений (t=0)

}

for(k=1;k<=Nt;k++)

{

solve(k*dt);

for(i=0;i<=n;i++)

{

u[i]=un[i];

if(i%10==0 && k%10==0) fprintf(fp,"%f\t%f\t%f\n",k*dt,i*h,u[i]); // печать в файл

}

}

fclose(fp);

}

void solve(double t) // основная процедура

{

double us0, us1;

un[0]=ux0(t); // u(0,t) - граничное условие

for(i=0; i<n; i++)

{

us1=u[i+1]; // начальное приближение

do { // цикл итерационного метода

us0=us1;

us1=us0-f(us0,un[i],u[i+1],u[i])/df(us0); // расчет следующего приближения

} while(abs(us1-us0)>eps); // выходим из цикла при достижении точности eps

un[i+1]=us1; // найденное значение

}

}

double g(double z)

{

return -log(exp(2.0*z)+1.0);

}

double dg(double z)

{

return -2.0*exp(2.0*z)/(1.0+exp(2.0*z));

}

double f(double z,double un0,double u1,double u0)

{

return (un0-u0+z-u1)/2.0/dt+(g(u1)-g(u0)+g(z)-g(un0))/2.0/h;

}

double df(double z)

{

return 1.0/2.0/dt+dg(z)/2.0/h;

}

double ut0(double z) // u(x,0) нач. усл.

{

return sin(pi*z);

}

double ux0(double z) // u(0,t) гр. усл.

{

return 0.0;

}