@XF 2016-11-20T16:08:02.000000Z 字数 2344 阅读 205

Exercise09

homewodk

exercise3.31

摘要

Study the tratrajectory of a small ball which can be considered as a point and will bounce at the boundary of a 'stadium' whose shape can be determined freely.

背景

Study the tratrajectory of a small ball which can be considered as a point and will bounce at the boundary of a 'stadium' whose shape can be determined freely.

正文

我们思考一个球在一个完美的台球桌上没有摩擦的问题，在碰撞之间的速度是恒定的，所以我们有

$\frac{\mathrm{d} x}{\mathrm{d} t}=v_{x}$

$\frac{\mathrm{d} y}{\mathrm{d} t}=v_{y}$
碰撞点确定后,接下来我们必须得到单位向量法在墙上的碰撞,

$\hat{n}$ 然后计算平行和垂直于墙的分量

$\hat{v_{i}}$ 我们有

$\vec{v}_{i,\perp }=(\vec{v}_{i}\cdot \hat{n})\hat{n}$

$\vec{v}_{i,\parallel }=\vec{v}_{i}-\vec{v}_{i,\perp }$
故从墙上反射会带的速度为

$\vec{v}_{i,f}=-\vec{v}_{i,\perp }$

$\vec{v}_{f,\parallel }=\vec{v}_{i,\parallel }$

结果

1.圆边界:

import matplotlib.pyplot as plt
import numpy as np
import math
class tabel:
    def__init__(self,x_0=0.2,y_0=0,vx_0=-0.2,vy_0=-0.5,dt_1=0.001,total_time=1000):
        self.x=[x_0]
        self.y=[y_0]
        self.vx=[vx_0]
        self.vy=[vy_0]
        self.dt=dt_1
        self.time=total_time
        self.t=[0]
    def calculate(self):
        for i in range(int(self.time/self.dt)):
            self.x.append(self.x[i]+self.vx[i]*self.dt)
            self.y.append(self.y[i]+self.vy[i]*self.dt)
            self.vx.append(self.vx[i])
            self.vy.append(self.vy[i])
            if ((self.x[i+1]**2+self.y[i+1]**2>1):
                X=(self.x[i]+self.x[i+1])/2
                Y=(self.y[i]+self.y[i+1])/2
                x1=X/(X**2+Y**2)**0.5
                y1=Y/(X**2+Y**2)**0.5
                self.vx[i+1]=(1-2*x1**2)*self.vx[i]-2*x1*y1*self.vy[i]
                self.vy[i+1]=(1-2*y1**2)*self.vy[i]-2*x1*y1*self.vx[i]
                if X**2+Y**2>1:
                    self.x[i+1]=X
                    self.y[i+1]=Y
                    continue
                else:
                    self.x[i]=X
                    self.y[i]=Y
                    continue
        x_1=[-1]
        y_1=[0]
        x_2=[-1]
        y_2=[0]
        dx=0.0001
    def bound(self):
        for k in range(20000):
            x_1.append(x_1[k]+dx)
            y_1.append((1-x_1[k+1]**2)**0.5)
            x_2.append(x_2[k]+dx)
            y_2.append(-(1-x_2[k+1]**2)**0.5)
    def show_results(self):
        plt.plot(self.x,self.y,'--',label='trajectory')
        plt.plot(x_1,y_1,'--',label='bound')
        plt.plot(x_2,y_2,'--',label='bound')
        plt.xlabel(u'x')
        plt.ylabel(u'y')
a=tabel()
a.calculate()
a.show_results()
b=tabel
b.bound()
b.show_results()
plt.show()

the results:for t=50,100,1000

the corresponding phase:

2.elliptical tabel:
The equation is:

$y^{2}+\frac{x^{2}}{4}=1$
the code which need to change:

if ((self.x[i+1]**2/4+self.y[I+1]**2)>1):
    X=(self.x[i]+self.x[i+1])/2
    Y=(self.y[i]+self.y[i+1])/2
    x1=(X/2)/(X**2/4+Y**2)**0.5
    x2=Y/(X**2/4+Y**2)**0.5
    self.vx[i+1]=(1-2*x1**2)*self.vx[i]-2*x1*y1*self.vy[i]
    self.vy[i+1]=(1-2*y1**2)*self.vy[i]-2*x1*y1*self.vx[i]