import math import pylab def drange start stop step start while stop a

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
import math
import pylab
def drange(start, stop, step):
l = []
r = start
while r <= stop:
l.append(r)
r += step
return l
A, B = 0.0 , 1.0
n = 5
h = (B - A) / n
D0, D1 = A + h, h
y = [[A, D0], [0, D1]]
def p(x): return 1
def q(x): return 1
def f(x): return x+1
def get_c1():
global n
return (B - y[0][n]) / y[1][n]
def get_solv_y_i(i): return y[0][i] + get_c1() * y[1][i]
x = drange(A, B, h)
def div(a, b):
return a / b
for i in range(1, n):
y[0].append(
div(
(h ** 2 * f(x[i]) - (1.0 - (h / 2) * p(x[i])) * y[0][i - 1] - (h ** 2 * q(x[i]) - 2) * y[0][i]),
1 + h / 2 * p(x[i])
)
)
y[1].append(
div(
-(1 - h / 2 * p(x[i])) * y[1][i - 1] - (h ** 2 * q(x[i]) - 2) * y[1][i],
1 + h / 2 * p(x[i])
)
)
pylab.plot(x, [get_solv_y_i(i) for i in range(n + 1)])
pylab.show()
for i in range(n):
print x[i], get_solv_y_i(i)