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

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

1 <= x <= 2 , y(1)= 0.5  ,  y(2)= 1/3

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


'''
/* ex7-4.py uses finite difference method to solve
 * 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.
 */
 '''
import math

def  p(x,y):
    return (y)
   
def q(x,y):
    return (0.0)

def r(x,y):
    return (math.pow(y,3))

def w(x):
    return (1.0 / (1+x) )


def tridg( n , a, b,  c,  d ):
    for i in range (2 , n+1):
        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 x[i] */
    d[n]=d[n]/b[n]

    for i in range (n-1 , 0 , -1):
        d[i]=(d[i]-c[i]*d[i+1]) / b[i]
    return

#++++++++++++++++++++++++++++++++++++++++
# n=99 aa=1.0 bb=2.0 alfa=0.5 bata=0.333333
#++++++++++++++++++++++++++++++++++++++++
n=99
aa=1.0
bb=2.0
alfa=0.5
bata=0.333333

a= [0.0 for i in range(0,n+1)]     #a [n] 矩陣
b= [0.0 for i in range(0,n+1)]     #a [n] 矩陣
c= [0.0 for i in range(0,n+1)]     #a [n] 矩陣
d= [0.0 for i in range(0,n+1)]     #a [n] 矩陣

y= [0.5 for i in range(0,n+1)]     #a [n] 矩陣


h=(bb-aa)/(n+1)

for k in range (1, 301):
    #x1=aa+h;
    x1=aa+h;
    b[1]=math.pow(h,2)* q ( (x1) , (y[1]) ) - 2.0
    c[1]=(1+(h/2.0)*p(  (x1) , (y[1]) ))
    d[1]=math.pow(h,2)*r( (x1) , (y[1]) ) - (1.0 - (h/2.0) * p(  (x1) , (y[1]))) * alfa

    for i in range (2 , n):
        x=aa+i*h;
        a[i]=1-(h/2.0)*p((x),(y[i]));
        b[i]=math.pow(h,2)*q((x),(y[i]))-2.0;
        c[i]=1+(h/2.0)*p((x),(y[i]));
        d[i]=math.pow(h,2)*r((x),(y[i]));

    xn=aa+n*h
    a[n]=1-(h/2.0)*p( (xn) , (y[n]) )
    b[n]=math.pow(h,2)*q((xn) , (y[n]))-2;
    d[n]=math.pow(h,2)*r((xn) , (y[n]))-(1+(h/2.0)*p( (xn) , (y[n]) ) )*bata;

    tridg(n,a,b,c,d)

    err=0.0
    for i in range (1 , n+1):
        err=err+abs(d[i]-y[i])
       
    if(err >0.001):
        for i in range (1, n+1):
            y[i]=d[i]
    else:
            break

#bound:
print("The iterations={%d}\n" %(k))
print(" x\t\ty(x)\t\t\tw(x)\t\t\t|y(x)-w(x)|\n")
print ("{%6.4f}\t\t{%10.7f}\t\t{%10.7f}\t\t{%10.7f}" %(aa,alfa, w(aa), abs(alfa-w(aa) ) ) )
for i in range (1 , n+1):
    x=aa+i*h;
    if(i%10==0):
        print ("{%6.4f}\t\t{%10.7f}\t\t{%10.7f}\t\t{%10.7f}" %(x,d[i],w(x),abs(d[i]-w(x) ) ))

print ("{%6.4f}\t\t{%10.7f}\t\t{%10.7f}\t\t{%10.7f}" %(bb,bata,w(bb),abs(bata-w(bb) ) ))

#for i in range (1,n+1) :      # //print the new matrix
 #   print( round(s1[i],4),"\t",end='')

輸出畫面
======== RESTART: F:\2018-09勤益科大數值分析\數值分析\PYTHON\EX7-4.py ===========
The iterations={5}

 x y(x) w(x) |y(x)-w(x)|

{1.0000} { 0.5000000} { 0.5000000} { 0.0000000}
{1.1000} { 0.4761904} { 0.4761905} { 0.0000000}
{1.2000} { 0.4545454} { 0.4545455} { 0.0000001}
{1.3000} { 0.4347825} { 0.4347826} { 0.0000001}
{1.4000} { 0.4166665} { 0.4166667} { 0.0000002}
{1.5000} { 0.3999998} { 0.4000000} { 0.0000002}
{1.6000} { 0.3846151} { 0.3846154} { 0.0000002}
{1.7000} { 0.3703701} { 0.3703704} { 0.0000003}
{1.8000} { 0.3571426} { 0.3571429} { 0.0000003}
{1.9000} { 0.3448273} { 0.3448276} { 0.0000003}
{2.0000} { 0.3333330} { 0.3333333} { 0.0000003}
>>>