2022全國大學生數學建模A題的思路與解法

import numpy as np

import matplotlib.pyplot as plt

from matplotlib.animation import FuncAnimation

import matplotlib

matplotlib.rcParams["font.sans-serif"] = ["SimHei"]

matplotlib.rcParams["axes.unicode_minus"] = False

def runge_kutta4(df, a, b, h, y0):

     num = len(y0)

     x = np.arange(a, b+h, h)

     w = np.zeros( (x.size, num) )

     w[0, :] = y0

     for i in range(x.size - 1):

         s0 = df(x[i], w[i, :],i*h)

         s1 = df(x[i] + h/2., w[i, :] + h * s0 / 2.,i*h)

         s2 = df(x[i] + h/2., w[i, :] + h * s1 / 2.,i*h)

         s3 = df(x[i+1], w[i, :] + h * s2,i*h)

         w[i+1,:] = w[i,:] + h * (s0 + 2*(s1+s2) + s3) / 6.

     return x, w

def df(x, variables,i):

     th1, th2, om1, om2 = variables

     A = np.zeros((2, 2))

     b = np.zeros(2)

     A[0, 0] = 2433+6201.535

     A[0, 1] = 2433

     A[1, 0] = 2433

     A[1, 1] = 2433

     b[0] = 6250*np.cos(1.4005*i)-10045*4*np.arctan(1)*th1*(th1<=0.999989731)-10045*4*np.arctan(1)*0.999989731*(th1>0.999989731)-656.3616*om1

     b[1] = -80000*th2-10000*om2

     dom1, dom2 = np.linalg.solve(A,b)

     return np.array([om1,om2,dom1,dom2])

a, b = 0.0,180.0

h = 0.01

th10 = 0.1

th20 = 0.1

om10 = 0.1

om20 = 0.1

y0 =跟單網gendan5.com np.array([th10, th20, om10, om20])

# 計算求解

t, w = runge_kutta4(df, a, b, h, y0)

th1 = w[:, 0]

th2 = w[:, 1]

om1 = w[:, 2]

om2 = w[:, 3]

plt.plot([h*i for i in range(len(th1))],th1,label="x1")

plt.plot([h*i for i in range(len(th2))],th2,label="x2")

#plt.plot([i*h for i in range(len(om1))],om1,label="x1'")

#plt.plot([i*h for i in range(len(om2))],om2,label="x2'")

plt.xlabel(" 時間 /s")

plt.ylabel(" 位移 /m",rotation=True)

#plt.ylabel(" 速度 /m*s^-1",rotation=True)

plt.legend()

plt.title("x1 x2 隨時間的變化量；求解步長 -m1"+str(h))

plt.savefig("x1 x2 隨時間的變化量 求解步長 -m1"+str(h)+".jpg")

plt.show()

#plt.title("x1' x2' 隨時間的變化量；求解步長 "+str(h))

#plt.savefig("x1' x2' 隨時間的變化量 求解步長 "+str(h)+".jpg")

#plt.show()

import xlsxwriter as xls

#th1 = np.array(th1[::20])

#th2 = np.array(th2[::20])

#om1 = np.array(om1[::20])

#om2 = np.array(om2[::20])

th1 = np.array(th1)

th2 = np.array(th2)

om1 = np.array(om1)

om2 = np.array(om2)

workbook = xls.Workbook("1-1-000-m1.xlsx")

worksheet = workbook.add_worksheet("Sheet1")

headings = ["x1","x2","v1","v2"]

worksheet.write_row("A1",headings)

worksheet.write_column("A2",th1)

worksheet.write_column("B2",om1)

worksheet.write_column("C2",th1+th2)

worksheet.write_column("D2",om1+om2)

workbook.close()