範例7-5 請用有限差法 (Finite difference Method) 解 非線性常微分方程式 y'' + y' ^ 2 + y = ln(x)

範例7-5 請用有限差法 (Finite difference Method) 解
非線性常微分方程式
y'' +  y' ^ 2  + y = ln(x)

1 <= x <= 2 , y(1)= 0.0  ,  y(2)= ln(2)=0.6931472

其真實解 W(x) = ln(x)
取 h=0.1  , n=99 , h=0.01 , n=99

/* ex7-5.c uses finite difference method to solve
 * nonlinear ordinary differential equation with boundary
 * conditions, y"+p(x)y'+q(x)y=r(x),a<=x<=b,
 * y(a)=alfa, y(b)=bata. After transfer ordinary
 * differential equation into system of linear algebra
 * equations, then call function tridg() to solve
 * tridiagonal equations.
 */
#include <stdio.h>
#include <math.h>
#define  p(x,y)  (y)
#define  q(x,y)  (1.0)
#define  r(x,y)  (log(x))
#define  w(x)    (log(x))
void tridg(int,double [],double [],double [],double []);
void main()
{
   int i,k,n;
   double a[100],b[100],c[100],d[100],y[100],dy[100],
  h,x1,x,xn,aa,bb,alfa,bata,err;
   scanf("n=%d aa=%lf bb=%lf alfa=%lf bata=%lf",
  &n,&aa,&bb,&alfa,&bata);
   h=(bb-aa)/(n+1);
   for(i=1;i<=n;i++)
      y[i]=0.0;
   
   for(k=1;k<=100;k++)
   {
        x1=aa+h;
        /* dy[1]=y1'=(y2-y0)/(2h) */
        dy[1]=(1.0/(2*h))*(y[2]-alfa);
        b[1]=pow(h,2)*q((x1),(y[1]))-2.0;
        c[1]=(1+(h/2.0)*p((x1),(dy[1])));
        d[1]=pow(h,2)*r((x1),(y[1]))-(1.0-(h/2.0)*
    p((x1),(dy[1])))*alfa;
        for(i=2;i<=n-1;i++)
        {
            x=aa+i*h;
        /* dy[i]=yi' */
        dy[i]=(1.0/(2*h))*(y[i+1]-y[i-1]);
            a[i]=1-(h/2.0)*p((x),(dy[i]));
            b[i]=pow(h,2)*q((x),(y[i]))-2.0;
            c[i]=1+(h/2.0)*p((x),(dy[i]));
            d[i]=pow(h,2)*r((x),(y[i]));
        }
        xn=aa+n*h;
        /* dy[n]=yn' */
        dy[n]=(1.0/(2*h))*(bata-y[n-1]);
        a[n]=1-(h/2.0)*p((xn),(dy[n]));
        b[n]=pow(h,2)*q((xn),(y[n]))-2.0;
        d[n]=pow(h,2)*r((xn),(y[n]))-(1+(h/2.0)*
        p((xn),(dy[n])))*bata;
        tridg(n,a,b,c,d);
        err=0.0;
        for(i=1;i<=n;i++)
        err=err+fabs(d[i]-y[i]);
        if(err >0.001)
        {
        for(i=1;i<=n;i++)
            y[i]=d[i];
        }
        else
        goto bound;
   }
   bound:
   printf("The iterations=%d\n",k);
   printf("x       y(x)       w(x)      |y(x)-w(x)|\n");
   printf("%5.3lf %10.7lf %10.7lf %10.7lf\n",
   aa,alfa,w(aa),fabs(alfa-w(aa)));
   for(i=1;i<=n;i++)
   {
      x=aa+i*h;
      if(i%10==0)
      printf("%5.3lf %10.7lf %10.7lf %10.7lf\n",
      x,d[i],w(x),fabs(d[i]-w(x)));
   }
   printf("%5.3lf %10.7lf %10.7lf %10.7lf\n",
    bb,bata,w(bb),fabs(bata-w(bb)));
   return;
}
void tridg(int n,double a[],double b[],double c[],double d[])
{
   int i;
   double r;
   for(i=2;i<=n;i++)
   {
      r=a[i]/b[i-1];
      b[i]=b[i]-r*c[i-1];
      d[i]=d[i]-r*d[i-1];
   }
   /* The answers are stored in d[i] */
   d[n]=d[n]/b[n];
   for(i=n-1;i>=1;i--)
      d[i]=(d[i]-c[i]*d[i+1])/b[i];
   return;
}


輸入資料
n=99 aa=1.0 bb=2.0 alfa=0.0 bata=0.6931472

輸出資料
The iterations=6
x       y(x)       w(x)      |y(x)-w(x)|
1.000  0.0000000  0.0000000  0.0000000
1.100  0.0953101  0.0953102  0.0000001
1.200  0.1823215  0.1823216  0.0000001
1.300  0.2623644  0.2623643  0.0000001
1.400  0.3364727  0.3364722  0.0000004
1.500  0.4054658  0.4054651  0.0000007
1.600  0.4700045  0.4700036  0.0000009
1.700  0.5306291  0.5306283  0.0000009
1.800  0.5877874  0.5877867  0.0000007
1.900  0.6418543  0.6418539  0.0000004
2.000  0.6931472  0.6931472  0.0000000