微分方程5_如何理解$e^i*pi

在之前的微积分学习中,通过实例了解到了 e t e^t et的性质

  1. e t e^t et的导数等于它本身,那么 d / d t ( e a t ) = a ∗ e a t d/dt(e^{at})=a*e^{at} d/dt(eat)=aeat

  2. 任何数的0次方等于1

在这里插入图片描述

假如函数 e t e^t et描述了数轴上的位置,它随着事件而变化y= e t e^t et,t代表时间,输出值y等于位置
在这里插入图片描述
那么位置将从1开始,因为t=0时y=1
在这里插入图片描述
因为 e t e^t et的导数等于本身,任何时候速度的值等于位置本身

在这里插入图片描述
通过位置与速度的关系,就可以只管的理解这个函数关系
在这里插入图片描述
当指数是常数倍的时候,速度就等于位置的常数倍
在这里插入图片描述
在这里插入图片描述
如果指数的倍数是负数,表示先翻转向量,再缩放
在这里插入图片描述
在这里插入图片描述

当倍数常数是代表-1开方的i,
在这里插入图片描述

在这里插入图片描述

在这里插入图片描述
旋转90°,已经超出了数轴,我们用复平面来表示

在这里插入图片描述

在这里插入图片描述

这个向量,就是 a + b ∗ i a+b*i a+bi

它的初始状态 e i ∗ 0.00 = 1.00 + 0.00. ∗ i e^{i*0.00}=1.00+0.00.*i ei0.00=1.00+0.00.i

在这里插入图片描述

1个单位时间后,速度为1个单位,走过了一个圆
在这里插入图片描述

π \pi π秒之后,到达了x轴负半轴,所以 e i π = − 1 + 0.00 i = − 1 e^{i\pi}=-1+0.00i=-1 eiπ=1+0.00i=1
在这里插入图片描述
τ \tau τ秒之后,( τ \tau τ=2 π \pi π),走了一个圆周的距离,
e i τ = 1 + 0.00 i = 1 e^{i\tau}=1+0.00i=1 eiτ=1+0.00i=1

在这里插入图片描述
由此可以看出,复指数可以看着做单位弧度
在这里插入图片描述

在这里插入图片描述

function dy = dynamic_system(t, y, wind_u , wind_v, current_u , current_v, spotld_lat) % 定义参数 rou_a=1.25;% 空气密度1.25kg/m3 Ca=0.7;% 风力系数 Sa=0.2;% 风作用在浮标上的面积,m3;pi*0.12^2 m = 6.35; % 物体质量 (kg) % 提取状态变量 u = y(3); v = y(4); [va, phi_a] = calculatewdir(wind_u , wind_v); % 计算风拖曳力 Fa_x = 0.5 * rou_a * Ca * Sa * (va * sin(phi_a) - u) * sqrt((va * sin(phi_a) - u)^2 + (va * cos(phi_a) - v)^2); Fa_y = 0.5 * rou_a * Ca * Sa * (va * cos(phi_a) - v) * sqrt((va * sin(phi_a) - u)^2 + (va * cos(phi_a) - v)^2); % 计算海流作用力需要的参数 rou_w=1020;% 海水密度1020kg/m3 Cw=1.6;% 流力系数,按照《船舶与海洋工程环境载荷》建议1.6 Sw=1;% 海流作用横截面积,m3;0.05*2*pi*0.075^2+pi*0.075^2+2*pi*0.152^2-2*pi*0.075^2 [vw, phi_w] = calculatewdir(current_u , current_v);%角度以y轴正向为零度,顺时针计算,范围0-2pi % 海流作用力 Fw_x = 0.5 * rou_w * Cw *Sw * (vw * sin(phi_w) - u) * sqrt((vw * sin(phi_w) - u)^2 + (vw * cos(phi_w) - v)^2); Fw_y = 0.5 * rou_w * Cw *Sw * (vw * cos(phi_w) - v) * sqrt((vw * sin(phi_w) - u)^2 + (vw * cos(phi_w) - v)^2); % 科氏力 omega = 7.929e-5;%度/小时 zeta = deg2rad(spotld_lat); f = 2 * omega * sin(zeta); Fc_x = - m * f * (u - vw * cos(phi_w)); Fc_y = - m * f * (v - vw * sin(phi_w)); % 计算状态变量的导数 dy = zeros(4, 1); dy(1) = u; dy(2) = v; dy(3) = Fa_x + Fw_x - Fc_x/ m; dy(4) = Fa_y + Fw_y - Fc_y/ m; end “主函数中调用” [tout1,yout1] = ode45(@(t, y) dynamic_system(t, y, wind_u , wind_v, current_u , current_v, spotld_lat), tspan, y0); 帮我改成python代码
最新发布
03-08
<think>嗯,用户让我把Matlab的代码转换成Python,特别是用ode45求解微分方程的部分。首先,我需要理解这段Matlab代码的功能。看起来这是一个动态系统模型,可能用于模拟浮标或其他海洋结构物在风力和海流作用下的运动。代码中的dynamic_system函数计算了状态变量的导数,包括位置和速度的变化,考虑了风拖曳力、海流作用力和科氏力。 接下来,我要考虑如何将Matlab代码转换为Python。Matlab的ode45在Python中通常用scipy的solve_ivp来实现,或者使用odeint,不过odeint现在推荐用solve_ivp替代。需要确认的是函数的参数传递方式,因为Matlab的ode45允许传递额外的参数,而Python的solve_ivp同样可以通过args参数传递。 首先,dynamic_system函数的参数在Matlab中是(t, y, wind_u, wind_v, current_u, current_v, spotld_lat),而Python的solve_ivp需要函数签名为func(t, y, *args),所以需要调整参数顺序,可能将额外的参数通过args传入。 然后,注意Matlab代码中的calculatewdir函数,这个函数可能需要自己实现,因为用户没有提供。这个函数应该将风速或海流的u和v分量转换为幅值和方向角。在Python中,可以用numpy的arctan2和hypot来计算角度和幅值。 接下来,处理变量提取部分。Matlab中的y(3)和y(4)对应Python中的y[2]和y[3],因为Python是0-based索引。需要调整所有索引。 然后,各个力的计算部分,比如Fa_x、Fa_y、Fw_x、Fw_y等,需要确保数学运算正确转换为Python的语法,比如平方根用np.sqrt,幂运算用**。 科氏力部分中的deg2rad可能需要用numpy的radians函数,或者math.radians,但需要注意spotld_lat是否为度数,这里Matlab的deg2rad将度数转换为弧度,Python中同样处理。 状态变量的导数dy在Matlab中是列向量,Python中可以返回列表,顺序对应dy[0], dy[1], dy[2], dy[3]。 主函数调用部分,ode45在Python中使用solve_ivp,需要设置t_span和y0,以及method='RK45'来对应ode45的算法。同时,args参数传递额外的参数,如wind_u, wind_v等。 需要注意的是,Matlab和Python在处理数组时的差异,比如Matlab的数组是列优先,而Python的numpy数组默认是行优先,但在这段代码中可能影响不大,因为主要是处理一维数组。 另外,需要确保所有变量和参数在Python中正确命名,比如rou_a改为rho_a,保持一致性,但用户可能希望保留原名,所以可能需要保留rou_a这样的变量名,或者根据习惯调整。 最后,测试代码是否正确,可能需要用户提供calculatewdir的实现,或者自己实现一个。例如,风速的u和v分量转换为幅值和方向角,可能需要计算va = sqrt(wind_u^2 + wind_v^2),方向角phi_a = arctan2(wind_u, wind_v),但需要确认Matlab中的方向定义,因为用户注释提到角度以y轴正向为零度,顺时针计算,范围0-2pi。这可能涉及到角度计算的方向和象限处理,需要仔细处理。 总结步骤: 1. 将dynamic_system函数转换为Python函数,注意参数顺序和索引。 2. 实现calculatewdir函数,将u和v分量转换为幅值和方向角。 3. 调整数学运算和numpy函数的使用。 4. 使用solve_ivp替换ode45,正确处理参数传递。 5. 测试代码,确保结果一致。</think>以下是转换后的Python代码实现,包含关键步骤说明: ```python import numpy as np from scipy.integrate import solve_ivp def calculatewdir(u, v): """将速度分量转换为幅值和方向角(弧度)""" magnitude = np.sqrt(u**2 + v**2) angle = np.arctan2(u, v) # 以y轴正向为零点,顺时针方向符合数学惯例 return magnitude, angle % (2*np.pi) # 确保角度在0-2&pi;范围内 def dynamic_system(t, y, wind_u, wind_v, current_u, current_v, spotld_lat): # 定义参数 rho_a = 1.25 # 空气密度(kg/m³) Ca = 0.7 # 风力系数 Sa = 0.2 # 风作用面积(m²) m = 6.35 # 质量(kg) rho_w = 1020 # 海水密度(kg/m³) Cw = 1.6 # 流力系数 Sw = 1.0 # 海流作用面积(m²) omega = 7.929e-5 # 地转角速度(rad/s) # 提取状态变量 [x, y, u, v] u = y[2] v = y[3] # 计算风参数 va, phi_a = calculatewdir(wind_u, wind_v) # 计算风拖曳力 rel_wind_u = va * np.sin(phi_a) - u rel_wind_v = va * np.cos(phi_a) - v Fa_x = 0.5 * rho_a * Ca * Sa * rel_wind_u * np.hypot(rel_wind_u, rel_wind_v) Fa_y = 0.5 * rho_a * Ca * Sa * rel_wind_v * np.hypot(rel_wind_u, rel_wind_v) # 计算海流参数 vw, phi_w = calculatewdir(current_u, current_v) # 计算海流作用力 rel_current_u = vw * np.sin(phi_w) - u rel_current_v = vw * np.cos(phi_w) - v Fw_x = 0.5 * rho_w * Cw * Sw * rel_current_u * np.hypot(rel_current_u, rel_current_v) Fw_y = 0.5 * rho_w * Cw * Sw * rel_current_v * np.hypot(rel_current_u, rel_current_v) # 计算科氏力 f = 2 * omega * np.sin(np.radians(spotld_lat)) # 转换为弧度 Fc_x = -m * f * (u - vw * np.cos(phi_w)) Fc_y = -m * f * (v - vw * np.sin(phi_w)) # 计算导数 dydt = np.zeros(4) dydt[0] = u # dx/dt dydt[1] = v # dy/dt dydt[2] = (Fa_x + Fw_x - Fc_x/m) / m # du/dt dydt[3] = (Fa_y + Fw_y - Fc_y/m) / m # dv/dt return dydt # 主函数调用示例 if __name__ == "__main__": # 示例参数 wind_u = 5.0 # 风速东向分量 wind_v = 3.0 # 风速北向分量 current_u = 0.5 # 海流东向分量 current_v = 0.3 # 海流北向分量 spotld_lat = 30.0 # 纬度(度) tspan = (0, 3600) # 时间范围(秒) y0 = [0, 0, 0, 0] # 初始状态 [x, y, u, v] # 求解ODE sol = solve_ivp( fun=dynamic_system, t_span=tspan, y0=y0, args=(wind_u, wind_v, current_u, current_v, spotld_lat), method='RK45', dense_output=True ) # 提取结果 t = sol.t y = sol.y.T ``` 主要修改说明: 1. **函数参数调整**- Python版本将`wind_u`, `wind_v`等参数通过`args`传递 - 状态变量`y`改为0-based索引(Matlab是1-based) 2. **数学计算改进**- 使用`np.hypot()`代替显式平方根计算,提高数值稳定性 - 使用`np.arctan2()`处理方向角计算,自动处理象限问题 3. **单位转换优化**- 使用`np.radians()`替代`deg2rad`进行角度转换 - 添加了角度归一化处理`angle % (2*np.pi)` 4. **数值积分器**- 使用`scipy.integrate.solve_ivp`代替`ode45` - 选择RK45方法(即ode45算法) 使用时需要注意: 1. 输入参数单位需保持一致(建议使用国际单位制) 2. 海流和风速分量需按东向(u)和北向(v)分量输入 3. 纬度参数`spotld_lat`以度为单位 4. 结果输出的时间步长由求解器自动确定,如需均匀间隔可后期插值 建议添加可视化代码观察结果: ```python import matplotlib.pyplot as plt plt.figure(figsize=(10,6)) plt.plot(y[:,0], y[:,1], label='Trajectory') plt.xlabel('X Position (m)') plt.ylabel('Y Position (m)') plt.title('Object Trajectory Simulation') plt.grid(True) plt.legend() plt.show() ```
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值