回答編集履歴
1
グラフ追加
answer
CHANGED
|
@@ -1,11 +1,10 @@
|
|
|
1
|
-
色々ややこしいので、間違いはご容赦ください。
|
|
2
1
|
x,yにリストとして座標を追加し、numpyの配列に変換すれば以下のソースコードで計算できると思います。
|
|
2
|
+
Ax=rとして解いた。a,b,c,dは三次関数の係数。
|
|
3
3
|
```python
|
|
4
4
|
import numpy as np
|
|
5
5
|
|
|
6
6
|
x = np.array([0, 1, 4, 5, 8])
|
|
7
7
|
y = np.array([0, 3, 3, 1, 2])
|
|
8
|
-
|
|
9
8
|
n = len(x)-1
|
|
10
9
|
h = x[1:]-x[:-1]
|
|
11
10
|
a = y
|
|
@@ -29,18 +28,23 @@
|
|
|
29
28
|
# Aの逆行列
|
|
30
29
|
Ainv = np.linalg.inv(A)
|
|
31
30
|
|
|
32
|
-
|
|
31
|
+
r = np.zeros(n+1)
|
|
33
|
-
b[:] = 0
|
|
34
|
-
|
|
32
|
+
r[1:-1] = (3/h)[1:]*(a[2:]-a[1:-1])-(3/h)[:-1]*(a[1:-1]-a[:-2])
|
|
35
33
|
|
|
36
|
-
c = np.dot(Ainv,
|
|
34
|
+
c = np.dot(Ainv, r)
|
|
37
35
|
|
|
36
|
+
d = (c[1:]-c[:-1])/(3*h)
|
|
37
|
+
b = (a[1:]-a[:-1])/h-h/3*(2*c[:-1]+c[1:])
|
|
38
38
|
print("x", x, "y", y)
|
|
39
39
|
print("h", h, "a", a)
|
|
40
40
|
print("A")
|
|
41
41
|
print(A)
|
|
42
|
+
print("r", r)
|
|
43
|
+
|
|
44
|
+
print("a", a)
|
|
42
45
|
print("b", b)
|
|
43
46
|
print("c", c)
|
|
47
|
+
print("d", d)
|
|
44
48
|
```
|
|
45
49
|
|
|
46
50
|
```text
|
|
@@ -52,6 +56,11 @@
|
|
|
52
56
|
[0. 3. 8. 1. 0.]
|
|
53
57
|
[0. 0. 1. 8. 3.]
|
|
54
58
|
[0. 0. 0. 0. 1.]]
|
|
55
|
-
|
|
59
|
+
r [ 0. -9. -6. 7. 0.]
|
|
60
|
+
a [0 3 3 1 2]
|
|
61
|
+
b [ 3.31018519 2.37962963 -1.96759259 -1.5462963 ]
|
|
56
62
|
c [ 0. -0.93055556 -0.51851852 0.93981481 0. ]
|
|
63
|
+
d [-0.31018519 0.04578189 0.48611111 -0.10442387]
|
|
57
|
-
```
|
|
64
|
+
```
|
|
65
|
+
|
|
66
|
+
|