題目回顧與寫(xiě)在前面

• 首先，我們隊(duì)在歷經(jīng)了千辛萬(wàn)苦之后，光榮得獲得了? 省三......
• 隊(duì)伍構(gòu)成 物理*2 + 計(jì)算機(jī)*1
• 隊(duì)伍分工? 計(jì)算機(jī)-->受力分析? 物理-->數(shù)值計(jì)算
• 總評(píng)：圖一樂(lè)，狠樂(lè)！物理系，計(jì)算機(jī)系嘛，不怎么看建模的啦！
  • 如果只是考慮力學(xué)問(wèn)題的話，我們分析得肯定還不太到位，但是這個(gè)題肯定很好分析
  • 所以，這個(gè)題主要在數(shù)值計(jì)算上

?思路

基本假設(shè)

1.海水是無(wú)粘及無(wú)旋的。

2.浮子在線性周期微幅波作用下會(huì)受到波浪激勵(lì)力（矩）、附加慣性力（矩）、興波阻尼力（矩）和靜水恢復(fù)力（矩）。

3.忽略中軸、底座、隔層及 PTO的質(zhì)量和各種摩擦。

4.初始浮子和振子平衡于靜水中。

問(wèn)題一的求解

?
問(wèn)題一浮子振子整體的受力分析	問(wèn)題一振子在浮子參考系下的受力分析

數(shù)值求解

4級(jí)顯式Runge-Kutta方法

求解過(guò)程的代碼：

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 = np.array([th10, th20, om10, om20])

# 計(jì)算求解
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("時(shí)間/s")
plt.ylabel("位移/m",rotation=True)
#plt.ylabel("速度/m*s^-1",rotation=True)
plt.legend()

plt.title("x1 x2 隨時(shí)間的變化量；求解步長(zhǎng) -m1"+str(h))
plt.savefig("x1 x2 隨時(shí)間的變化量 求解步長(zhǎng) -m1"+str(h)+".jpg")
plt.show()


#plt.title("x1' x2' 隨時(shí)間的變化量；求解步長(zhǎng) "+str(h))
#plt.savefig("x1' x2' 隨時(shí)間的變化量 求解步長(zhǎng) "+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()

問(wèn)題二

變步長(zhǎng)搜索法

矩形法數(shù)值積分

問(wèn)題二解法

60秒-80秒	80秒-100秒	100秒-120秒

?

?

問(wèn)題三

?

數(shù)值解（類(lèi)似一，二的解法）

?

?問(wèn)題四

?我們的優(yōu)化模型：

文章來(lái)源地址http://www.zghlxwxcb.cn/news/detail-455358.html

小結(jié)

• 這道題受力分析的成分要更大一點(diǎn)哈
• 分析完了就數(shù)值計(jì)算
• 受力分析不完整，我們就不放出來(lái)了，以后再完善吧。

到了這里，關(guān)于2022全國(guó)大學(xué)生數(shù)學(xué)建模A題的思路與解法的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！