# Position Based Dynamic 小练习 PS: 请问有没有同学用Taichi 实现过Projective Dynamics呀 3 Likes

I got a small program that connects springs to form something like it:

``````import taichi as ti
ti.init()

N = 8
x = ti.Vector.field(2, float, N, needs_grad=True)  # position
v = ti.Vector.field(2, float, N)  # velocity
U = ti.field(float, (), needs_grad=True)  # square error
c = ti.Vector.field(2, float, 1)

@ti.kernel
def compute_U():
for i in ti.ndrange(N-1):
r = x[i] - x[i+1]
U[None] += k * abs(r.norm() - dx)

@ti.kernel
for i in x:
v[i] += dt * (-x.grad[i] - friction * v[i] - g)  # dv/dt = -dU/dx, with added friction
for i in x:
x[i] += dt * v[i]  # dx/dt = v
x = x0_pos

def substep():
with ti.Tape(U):
# every kernel invocation within this indent scope
# will also be accounted into the partial derivative of U
# with corresponding input variables like x.
compute_U()  # will also computes dU/dx and save in x.grad

@ti.kernel
def init():
for i in range(N):
x[i] = x0_pos + i * c

# hooke's constant
k = 0.19
dx = 0.1
dx_sq = dx*dx

# gravity acceleration, vertical
g = [0., 0.01]

friction = [0.015, 0.08]
c = [0.1, -0.008]
x0_pos = [0.4, 0.9]
dt = 1e-3

it = 0
debug = 2

init()

gui = ti.GUI('Hanging-Springs')
while gui.running:
for i in range(50):
substep()
if debug == 1:
print('U = ', U[None])
elif debug == 2:
max_dx_sq = 0.
q_norm = 0.
q = [0., 0.]
p = x.to_numpy()
for i in range(N-1):
q = p[i] - p[i+1]
q_norm = q**2 + q**2
if q_norm > max_dx_sq:
max_dx_sq = q_norm
print("max diff = ", max_dx_sq - dx_sq)
it += 1