10.2 陰解法

// heat1d-i-freefem.edp
// 実行: FreeFem++ heat1d-i-freefem.edp
// N, x が大域的な識別子を隠すと警告が出る (つまり名前が衝突する)。

int i, N, n,n2, nMax,nfunc,nm1;
real pi=4.0*atan(1.0);
real a = 0.0, b = 1.0;
real lambda,Tmax,theta,c1,c2,c3,c4,dt,t;
int skip;

cout << "input nfunc= ";
cin >> nfunc;
cout << "input theta= ";
cin >> theta;
cout << "input N= ";
cin >> N;
cout << "input lambda= ";
cin >> lambda;
cout << "input Tmax= ";
cin >> Tmax;
cout << "input dt= ";
cin >> dt;


real h = (b-a) / N;  
real tau = lambda * h^2;
real[int] x(N+1);
real[int] u(N+1);
real[int] AL(N),AD(N),AU(N),ff(N);
real[int] al(N),ad(N),au(N);
real[int] al2(N),ad2(N),au2(N),b2(N);

skip = rint(dt / tau);

cout << "h= " << h << endl;
cout << "lambda= " << lambda << endl;
cout << "tau= " << tau << endl;

c1 = - theta*lambda;
c2 = 1.0 + 2.0*theta*lambda;
c3 = 1.0 - 2.0*(1.0-theta)*lambda;
c4 = (1.0 -theta)*lambda;

for (i = 1; i <N;i++){
  AL[i] = AU[i] = -theta*lambda;
  AD[i] = 1 +2*theta*lambda;
}

func real f(real x) {
  if (nfunc ==1)
    if (x < 0.5)
      return x;
    else
      return 1 - x;
  else if (nfunc == 2)
    return 1.0;
  else if (nfunc == 3)
    return sin(pi*x);
  else if (nfunc == 4)
    if (x <=0.1)
      return 5*x;
    else if (x <= 0.3)
      return -2*x + 0.7;
    else if (x <= 0.5)
      return 4.5*x - 1.25;
    else if (x <= 0.7)
      return -x + 1.5;
    else if (x <= 0.9)
      return x + 0.1;
    else
      return -10*x + 10.0;
  else if (nfunc == 5)
    if (x <= 0.2)
      return -x + 0.8;
    else if (x <= 0.3)
      return 4*x - 0.2;
    else if (x <= 0.4)
      return -4*x + 2.2;
    else if (x <= 0.7)
      return -x + 1.0;
    else if (x <= 0.8)
      return 4*x - 2.6;
    else 
      return -x +1.5;
  else
    return 0;
}

for (i = 0; i <= N; i++) {
  x[i] = a + i * h;
  u[i] = f(x[i]);
}
plot([x,u], bb=[[-0.1,-0.1],[1.1,1.1]],aspectratio=true);

n2 = N-1;
for (i = 1; i <N;i++){
  al[i] = AL[i]+1;
  ad[i] = AD[i]+1;
  au[i] = AU[i]+1;
}
for (i = 1; i <N;i++){
  al2[i] = AL[i]+1;
  ad2[i] = AD[i]+1;
  au2[i] = AU[i]+1; 
}
//係数行列の LU 分解
nm1 = n2 -1;
for (i=0; i < nm1; i++){
    al[i+1] /= au[i];
    ad[i+1] -= au[i]*al[i+1];
}

nMax = rint(Tmax / tau);
cout << "nMax =" << nMax << endl;
for (n = 1; n <= nMax; n++) {
//連立1次方程式を作って解く
  for (i = 1; i < N; i++){
    ff[i] = c3*u[i] +c4*(u[i-1] +u[i+1]);
  }
// LU 分解済みの三重対角行列を係数に持つ三項方程式を解く
  for(i=0; i < nm1;i++){
    b2[i] = ff[i]+1;
    b2[i+1] -= b2[i] * al2[i+1];
  }
  b2[nm1]/=ad2[nm1];
  for(i=n2-2; i >= 0 ;i-1){
    b2[i] = (b2[i] - au2[i]*b2[i+1])/ad2[i];
  }
  for(i=1;i<N;i++){
    u[i] =b2[i];
  }
  u[0]=0.0;
  u[N]=0.0;
  t=n*tau;
  if(n%skip ==0){
    plot([x, u], bb=[[-0.1,-0.1],[1.1,1.1]],aspectratio=true);
  }
}