numerical-analysis-python/ex8.1.py at main · akshaycodeforge/numerical-analysis-python · GitHub

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
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Dec 14 19:06:20 2019

@author: astro
"""

from math import sin
import numpy as np
import matplotlib.pyplot as plt

# for legend
fig = plt.figure()
ax = plt.subplot(111)

#Define function
def f(x,t):
    return -x**3+sin(t)

a=0.0    #start of interval
b=100.0   #End of interval
N=1000    #Number of steps
h=(b-a)/N  #size of a single step


################################### Euler method


euler =np.arange(a,b,h)
eul=[]
x=0.0 #initial condition
for t in euler:
    eul.append(x)
    x+=h*f(x,t)
ax.plot(euler,eul,'k+',label='Euler')

######################################rk2 method


srk =np.arange(a,b,h)
mid = []
x=0.0 #initial condition
for t in srk:
    mid.append(x)
    k1 = h*f(x,t)
    k2 = h*f(x+0.5*k1,t+0.5*h)
    x += k2
ax.plot(srk,mid,'g*',label='Midpoint')

######################################## Rk4 method

tpoints =np.arange(a,b,h)
xpoints = []
x = 0.0 #initial condition
for t in tpoints:
    xpoints.append(x)
    k1 = h*f(x,t)
    k2 = h*f(x+0.5*k1,t+0.5*h)
    k3 = h*f(x+0.5*k2,t+0.5*h)
    k4 = h*f(x+k3,t+h)
    x += (k1+2*k2+2*k3+k4)/6
ax.plot(tpoints,xpoints,'r-',label='RK4')

#### for the graph

plt.xlabel("t")
plt.ylabel("x(t)")
ax.legend()
ax.legend(loc='upper right')
plt.show()