例題7-3 常微分方程式邊界 y'' + x^2 y' + xy = 2 + 3x^3 + (1+x+x^2) exp(x)

例題7-3 常微分方程式邊界

y'' + x^2 y' + xy = 2 + 3x^3 + (1+x+x^2) exp(x)

0<= x <= 1
y(0) = 1
y(1)= 3.7182818

h=0.1 或 h=0.01

'''
/* ex7-3.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):
    return (math.pow(x,2))
    
def q(x):
    return (x)

def r(x):
    return (3*math.pow(x,3)+math.exp(x)*(1.0+x+math.pow(x,2))+2.0)

def w(x):
    return (math.pow(x,2)+math.exp(x))


def tridg( n , a, b,  c,  d , s):
    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] */
    s[n]=d[n]/b[n]

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

    return s;

#++++++++++++++++++++++++++++++++++++++++
n=99
aa=0.0
bb=1.0
alfa=1.0
bata=3.7182818
a= [0.0 for i in range(n+1)]     #a [n] 矩陣
b= [0.0 for i in range(n+1)]     #a [n] 矩陣
c= [0.0 for i in range(n+1)]     #a [n] 矩陣
d= [0.0 for i in range(n+1)]     #a [n] 矩陣
s= [0.0 for i in range(n+1)]     #a [n] 矩陣
s1= [0.0 for i in range(n+1)]     #a [n] 矩陣

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


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

s1=tridg(n,a,b,c,d,s)

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,s1[i],w(x),abs(s1[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) ) ))






輸出畫面
======== RESTART: F:\2018-09勤益科大數值分析\數值分析\PYTHON\EX7-3.py ==========
 x y(x) w(x) |y(x)-w(x)|

{0.0000} { 1.0000000} { 1.0000000} { 0.0000000}
{0.1000} { 1.1151718} { 1.1151709} { 0.0000008}
{0.2000} { 1.2614043} { 1.2614028} { 0.0000016}
{0.3000} { 1.4398610} { 1.4398588} { 0.0000022}
{0.4000} { 1.6518274} { 1.6518247} { 0.0000027}
{0.5000} { 1.8987242} { 1.8987213} { 0.0000030}
{0.6000} { 2.1821218} { 2.1821188} { 0.0000030}
{0.7000} { 2.5037555} { 2.5037527} { 0.0000028}
{0.8000} { 2.8655432} { 2.8655409} { 0.0000023}
{0.9000} { 3.2696045} { 3.2696031} { 0.0000014}
{1.0000} { 3.7182818} { 3.7182818} { 0.0000000}
>>>