python物理建模初学者指南第四章,Python物理建模初学者指南(part II)

Scripts templete

#!/usr/bin/env python3

# -*- coding: utf-8 -*-

"""

Created on %(date)s

@author: %(ZYJ)s (Python 3.7)

@email: %(lianlizyj@mail.dlut.edu.cn)s

Description:

"""

import numpy as np

import matplotlib.pyplot as plt

#s=str 3 # SyntaxError

#prnt(s) # NameError

"""

for i in np.arange(-1,8):

print(i,np.log(i))

RuntimeWarning: divid by zero encontered in log

"""

"""

same as error function in matlab

for i in np.arange(-1,8):

assert (i>0),"I don't know how to take the log of {}!".format(i) # AssertionEeeor

print(i,np.log(i))

"""

# first_script.py

# Jesse M. Kinder -- 2014 Nov 11

"""

Calculate how long an object is in the air when thrown from a specied height

with a range of initial speeds assuming constant acceleration due to gravity:

o.5*g*t**2-v0*t-y0=0

"""

# %% Initialize variables.

ini_speed=0.0 # v0 = initial vertical speed of ball [m/s]

impact_time= 0.0 # t = time of impact [s] (computered in loop)

# %% Initialize parameters.

g = 9.8066 # Gravitational acceleration [m/s^2]

ini_height = 2.0 # y0 = height ball is thrown from [m]

speed_increament = 5.0 # Speed increment for each iteration [m/s]

cutoff_time = 10.0 # Stop computing after impact time exceeds cutoff

# %% Calculate and display impact time. Incerement initial speed each step.

# Repeat until impact time exceeds cutoff

while impact_time < cutoff_time:

# Use quadratic equation to solve kinematic equation for impact time:

impact_time = (np.sqrt(ini_speed**2+2*g*ini_height)\

+ini_speed)/g

print("Speed = {} m/s; time= {:.1f} s".format(ini_speed,impact_time))

ini_speed+=speed_increament

print("Calculation complete.")

# %% comments for variables define and unit

# Variables:

# ------------------------------------------

# length = length of the microtubule [um]

# velocity = velocity of motor [um/s]

# rate_constant = rate constant [1/days]

或有行为

#!/usr/bin/env python3

# -*- coding: utf-8 -*-

"""

Created on Wed Aug 26 20:35:33 2020

@author: Yanjie Zhang (Python 3.7)

@email: lianlizyj@mail.dlut.edu.cn

Description:

Purpose:

This script illustrates branching.

"""

import numpy as np

import matplotlib.pyplot as plt

for trial in range(5):

userInput=input('Pick a number: ')

number=float(userInput)

if number<0:

print('Square root is not real')

else:

print('Square root of {} is {:.4f}'.format(number,np.sqrt(number)))

userAgain=input('Another [y/n]? ')

if userAgain!='y':

break

# %% Check whether the loop terminated normally

if trial==4:

print('Sorry, only 5 per customer.')

elif userAgain=='n':

print("Bye!")

else:

print('Sorry, i didn\'t understand that.')

# %% nested loops

rows=3

columns=4

a=np.zeros((rows,columns))

for i in range(rows):

for j in range(columns):

a[i,j]=i**2+j**3

print("New a is : ", a)

数据输入、结果输出

#!/usr/bin/env python3

# -*- coding: utf-8 -*-

"""

Created on Wed Aug 26 21:36:28 2020

@author: Yanjie Zhang (Python 3.7)

@email: lianlizyj@mail.dlut.edu.cn

Description:

Purpose:

import data to python from *,csv(separator is ',' )

"""

import numpy as np

import matplotlib.pyplot as plt

path='HIVseries.csv'

data_set=np.loadtxt(path,delimiter=',')

# %% other data type process (open method)

my_file=open(path)

temp_data=[]

for line in my_file:

print(line)

x,y=line.strip().split(',')

temp_data=[float(x),float(y)]

my_file.close()

data_set=np.array(temp_data)

# %% with open method

with open(path) as my_file:

for line in my_file:

print(line)

x,y=line.strip().split(',')

temp_data=[float(x),float(y)]

data_set=np.array(temp_data)

# %% get data_set from web

# import urllib

# web_file=urllib.request.urlopen('http://www.physics.upenn.edu/biophys/'+\

# 'PMLS/Datasets/01HIVseries/HIVseries.csv')

# print(" read data from web")

# print(web_file)

# data_set=np.loadtxt(web_file)

# url='http://www.baidu.com'

# content=urllib.request.urlopen(url).read()

# print(content)

# %% save data in files

x=np.linspace(0,1,1001)

y=3*np.sin(x)**3-np.sin(x)

np.savetxt('x_value.dat',x) # readable

np.savetxt('y_value.dat',y)

np.save('x_values',x) # unreadable

np.save('y_values',y)

np.savez('xy_values',x_vals=x,y_vals=y) # unreadable

print('save data in file sucessfully!')

# %% load the data

y2=np.loadtxt('y_value.dat')

print('read data by loadtxt success!',y2)

x2=np.load('x_values.npy')

print("read data by load success!",x2)

w=np.load('xy_values.npz')

print('read data of *.npz by load success!',w)

print("read data of *.npz by key, key() is ",w.keys())

print('w[\'x_vals\'] is',w['x_vals'])

# %% print into file by open function

with open('power.txt','w') as my_file:

print(" N \t\t2**N\t\t3**N\n") # print labels for column

print("---\t\t----\t\t----\n") # print separator.

my_file.write(" N \t\t2**N\t\t3**N\n") # write labels to file

my_file.write("---\t\t----\t\t----\n") # write separator to file

# %% Loop over integerfrom 0 to 10 and print/write results.

for N in range(11):

print("{:d}\t\t{:d}\t\t{:d}\n".format(N,pow(2,N),pow(3, N)))

my_file.write("{:d}\t\t{:d}\t\t{:d}\n".format(N,pow(2,N),pow(3, N)))

re_read=np.loadtxt('power.txt',skiprows=2)

print(re_read)

print("How to read power.txt, this is a question!")

# read by specifing the dtype

re_read=np.loadtxt('power.txt',\

dtype={'names':('col1','col2','col3')\

,'formats':('float','float','float')},\

skiprows=2)

print(list(re_read))

print("Sorry, i can't read complex format txt file now.")

print("Now, i can't read txt by formats type. Yeah!")

# re_read=np.fromregex()

print("I think i don't need to use fromregex to read complex data type, now.")

数据可视化

#!/usr/bin/env python3

# -*- coding: utf-8 -*-

"""

Created on Wed Aug 26 23:52:26 2020

@author: Yanjie Zhang (Python 3.7)

@email: lianlizyj@mail.dlut.edu.cn

Description:

Purpose:

Learning to visualize the data

"""

import numpy as np

import matplotlib.pyplot as plt

from matplotlib.ticker import (MultipleLocator,FormatStrFormatter,\

AutoMinorLocator)

num_points=50

x_min,x_max=0,4

x_values=np.linspace(x_min, x_max,num_points)

y_values=x_values**2

# check for the length of x and y

assert len(x_values)==len(y_values),\

"Length-mismtch: {:d} versus {:d}".format(len(x_values),len(y_values))

plt.plot(x_values, y_values,'r--d')

# plt.scatter(x_values, y_values)

# plt.semilogy(x_values, y_values,'r-o')

# plt.loglog(x_values, y_values,'r-o')

plt.xlim(1,3) # xlim([1,3]) in matlab

plt.axis('auto') # axis equal/square/tight/auto in matlab

ax=plt.gca()

ax.set_title("My first plot in python",size=24, weight='bold',family='Times New Roman',fontsize=30,fontweight=1.5)

ax.set_xlabel("x [$\mu$m]",family='Times New Roman', weight='bold',fontsize=15,fontweight=1.5)

ax.set_ylabel('y',family='Times New Roman', weight='bold',fontsize=15,fontweight=1.5)

lines=ax.get_lines()

plt.setp(lines[0],linestyle='--',linewidth=3,color='c') # set_property

# %% use "label" keyword to set labels when plotting

plt.plot(x_values,y_values,label="population 1")

plt.plot(x_values,x_values**3,label= "population 2")

# plt.plot(x_values,y_values)

# plt.plot(x_values,x_values**3)

plt.legend()

# use line objects to set labels after plotting.

ax=plt.gca()

lines=ax.get_lines()

lines[0].set_label("Infected population")

lines[1].set_label("Cured population")

ax.legend()

# =============================================================================

# Make a plot with major ticks that are multiples of 20 and minor ticks that

# are multiples of 5. Label mahor ticks with '%d' formatting but don't label

# minorticks.

# =============================================================================

# ax.xaxis.grid(True,which='minor')

# ax.yaxis.grid(True,which='minor')

ax.xaxis.set_major_locator(MultipleLocator(0.5))

ax.xaxis.set_major_formatter(FormatStrFormatter('%.2f'))

ax.yaxis.set_major_locator(MultipleLocator(10))

ax.yaxis.set_major_formatter(FormatStrFormatter('%d'))

# =============================================================================

# For the minor ticks, use no labels; default NullFormatter

# =============================================================================

ax.xaxis.set_minor_locator(MultipleLocator(0.1))

# ax.yaxis.set_minor_locator(MultipleLocator(2))

ax.yaxis.set_minor_locator(AutoMinorLocator())

# change the orientation of ticks (in/out/inout)

plt.rcParams['xtick.direction']='in'

plt.rcParams['ytick.direction']='in'

ax.tick_params(which='both',width=2)

ax.tick_params(which='major', length=7,color='b')

ax.tick_params(which='minor', length=4,color='r')

# change the tick label property

plt.xticks(family='Times New Roman', weight='normal',fontsize=12)

plt.yticks(family='Times New Roman', weight='normal',fontsize=12)

# You can specify a rotation for tick labels in degrees or with keywords

plt.yticks(rotation='vertical')

# =============================================================================

# # Use Axes object to set labels after plotting (cmd+4)

# =============================================================================

ax.legend(("Healthy","Recovered"),loc=0) # upper/lower right/left/center

plt.show()

# %% plot with errorbar

plt.figure()

x_errors=np.random.random(len(x_values))

y_errors=np.random.random(len(x_values))

# plt.errorbar(x_values,y_values,yerr=y_errors,xerr=x_errors)

plt.errorbar(x_values,y_values,yerr=y_errors,fmt='o',ecolor='r',color='b',elinewidth=2,capsize=4)

plt.show()

# %% plot 3D graphy

from mpl_toolkits.mplot3d import Axes3D

fig=plt.figure() # create a new figure

ax=Axes3D(fig) # create 3D plotterattached to figure

t=np.linspace(0,5*np.pi,501) # define parameter for parameteric plot

ax.plot(np.cos(t),np.sin(t),t) # draw 3D plot

ax.view_init(elev=30,azim=30) # change the view

# %% multiple plot

x=np.linspace(0,1,5)

y1=np.exp(x)

y2=x**2

fig=plt.figure()

plt.plot(x,y1,'r--^',x,y2,'k-o')

# =============================================================================

# test: 3B

# =============================================================================

num_curves=3

x=np.linspace(0,1,501)

y=np.zeros((x.size,num_curves))

for n in range(num_curves):

y[:,n]=np.sin((n+1)*x*2*np.pi)

plt.figure()

plt.plot(x,y)

plt.legend(("2$\pi$","4$\pi$","6$\pi$"),loc='best')

# =============================================================================

# cla & hold

# =============================================================================

# fig=plt.figure()

t=np.linspace(0,5*np.pi,501)

plt.plot(t,np.sin(t)**2)

plt.plot(t,np.sin(t)**3)

ax=plt.gca()

plt.cla()

# ax.hold(True)

print("Hold attribute has been removed from matplotlib 3")

plt.plot(t,np.sin(t)**4)

print("You can close figure by plt.close() or plt.close('name') or plt.close('all')")

# %% subplot

from numpy.random import random

t=np.linspace(0,1,101)

plt.figure()

plt.subplot(2,2,1);plt.hist(random(20))

plt.subplot(2,2,2);plt.plot(t,t**2,t,t**3-t)

plt.subplot(2,2,3);plt.plot(random(20),random(20),'r*')

plt.subplot(2,2,4);plt.plot(t*np.cos(10*t),t*np.sin(10*t))

plt.tight_layout()

# fig=plt.gcf()

plt.savefig("greatest_figure_ever.pdf")

function

# -*- coding: utf-8 -*-

"""

Created on Sat Aug 29 15:10:44 2020

@author: Yanjie Zhang

@email: lianlizyj@mail.dlut.edu.cn

Description:

Purpose:

"""

import numpy as np

import matplotlib as plt

from matplotlib.ticker import (MultipleLocator,FormatStrFormatter,\

AutoMinorLocator)

def taxicab(pointA,pointB):

"""

Taxicab metric for computing distance between points A and B

pointA=(x1,y1)

pointB=(x2,y2)

Return |x2-x1|+|y2-y1|, Distances are measured in city blocks.

"""

interval=abs(pointB[0]-pointA[0])+abs(pointB[1]-pointA[1])

return interval

# %%

fare_rate=0.4

start=(1,2)

stop = (4,5)

a=taxicab(start, stop) *fare_rate

print(a)