题目七解答说明

y₁' = y₂y₃
y₂' = -y₁y₃  
y₃' = 0.51y₁y₂

dS/dt = -β(SI/N)
dI/dt = β(SI/N) - γI
dR/dt = γI

using DifferentialEquations
using Printf

println("题目七解答 - 使用ode45函数")
println("="^50)

# 问题1：刚体运动微分方程组
println("问题1：刚体运动微分方程组")
println("方程组：")
println("y₁' = y₂y₃")
println("y₂' = -y₁y₃")
println("y₃' = 0.51y₁y₂")
println("初始条件：y₁(0) = 0, y₂(0) = 1, y₃(0) = 1")
println()

# 刚体运动微分方程组 - ode45风格
function rigid_body_system!(du, u, p, t)
    y1, y2, y3 = u
    
    du[1] = y2 * y3      # y₁' = y₂y₃
    du[2] = -y1 * y3     # y₂' = -y₁y₃
    du[3] = 0.51 * y1 * y2  # y₃' = 0.51y₁y₂
end

# 初始条件
u0 = [0.0, 1.0, 1.0]
tspan = (0.0, 12.0)  # 修改时间范围为0到12，与图片一致

# 创建ODE问题
prob = ODEProblem(rigid_body_system!, u0, tspan)

# 使用ode45求解器 (DP5就是Julia版本的ode45)
sol = solve(prob, DP5(), reltol=1e-8, abstol=1e-8)

# 提取结果
t_values = sol.t
y1_values = [sol[i][1] for i in 1:length(sol)]
y2_values = [sol[i][2] for i in 1:length(sol)]
y3_values = [sol[i][3] for i in 1:length(sol)]

println("刚体运动方程组求解结果：")
println("时间范围：0 到 12")
println("最终值：")
@printf("y₁(12) = %.6f\n", y1_values[end])
@printf("y₂(12) = %.6f\n", y2_values[end])
@printf("y₃(12) = %.6f\n", y3_values[end])
println()

# 验证守恒量
println("验证守恒量：")
conserved_initial = y1_values[1]^2 + y2_values[1]^2 + y3_values[1]^2
conserved_final = y1_values[end]^2 + y2_values[end]^2 + y3_values[end]^2
@printf("初始值：y₁² + y₂² + y₃² = %.6f\n", conserved_initial)
@printf("最终值：y₁² + y₂² + y₃² = %.6f\n", conserved_final)
@printf("变化量：%.6f\n", abs(conserved_final - conserved_initial))
println()

println("刚体运动方程数值解计算完成！")
println("数据已保存到文件，可在Syslab中绘制图形。")
println()

println("="^50)

# 问题2：SIR传染病模型
println("问题2：SIR传染病模型")
println("模型参数：")
println("总人口：N = 10000")
println("初始感染者：I(0) = 10")
println("感染率：β = 0.4 day⁻¹")
println("恢复率：γ = 0.1 day⁻¹")
println("时间范围：0 到 100 天")
println()

# SIR模型微分方程组 - ode45风格
function sir_system!(du, u, p, t)
    S, I, R = u
    N = 10000.0
    β = 0.4
    γ = 0.1
    
    du[1] = -β * S * I / N  # dS/dt = -β(SI/N)
    du[2] = β * S * I / N - γ * I  # dI/dt = β(SI/N) - γI
    du[3] = γ * I  # dR/dt = γI
end

# 初始条件
N = 10000.0
I0 = 10.0
S0 = N - I0
R0 = 0.0

u0_sir = [S0, I0, R0]
tspan_sir = (0.0, 100.0)

# 创建ODE问题
prob_sir = ODEProblem(sir_system!, u0_sir, tspan_sir)

# 使用ode45求解器
sol_sir = solve(prob_sir, DP5(), reltol=1e-8, abstol=1e-8)

# 提取结果
t_sir_values = sol_sir.t
S_values = [sol_sir[i][1] for i in 1:length(sol_sir)]
I_values = [sol_sir[i][2] for i in 1:length(sol_sir)]
R_values = [sol_sir[i][3] for i in 1:length(sol_sir)]

println("SIR模型求解结果：")
println("初始状态：")
@printf("易感者 S(0) = %.0f\n", S_values[1])
@printf("感染者 I(0) = %.0f\n", I_values[1])
@printf("康复者 R(0) = %.0f\n", R_values[1])
println()

println("最终状态（100天后）：")
@printf("易感者 S(100) = %.0f\n", S_values[end])
@printf("感染者 I(100) = %.0f\n", I_values[end])
@printf("康复者 R(100) = %.0f\n", R_values[end])
println()

# 找到感染高峰
max_infected, max_index = findmax(I_values)
max_time = t_sir_values[max_index]
@printf("感染高峰：第 %.1f 天，感染者数量 %.0f\n", max_time, max_infected)
println()

# 验证总人口守恒
total_initial = S_values[1] + I_values[1] + R_values[1]
total_final = S_values[end] + I_values[end] + R_values[end]
@printf("总人口守恒验证：")
@printf("初始总人口：%.0f\n", total_initial)
@printf("最终总人口：%.0f\n", total_final)
@printf("误差：%.6f\n", abs(total_final - total_initial))
println()

# 计算基本再生数 R₀
β = 0.4
γ = 0.1
R0_basic = β / γ
@printf("基本再生数 R₀ = β/γ = %.1f\n", R0_basic)

# 计算最终易感者比例
final_susceptible_ratio = S_values[end] / 10000
@printf("最终易感者比例：%.4f\n", final_susceptible_ratio)

# 计算总感染人数
total_infected = 10000 - S_values[end]
@printf("总感染人数：%.0f\n", total_infected)
@printf("总感染比例：%.4f\n", total_infected / 10000)

println()
println("="^50)
println("解答完成！")

# 输出一些关键时间点的数据
println()
println("关键时间点数据：")
println("时间(天) | 易感者 | 感染者 | 康复者")
println("-"^40)

# 找到关键时间点的数据
key_times = [1, 10, 20, 30, 50, 70, 100]
for time in key_times
    # 找到最接近的时间点
    idx = argmin(abs.(t_sir_values .- time))
    @printf("%7.0f | %6.0f | %6.0f | %6.0f\n", 
            t_sir_values[idx], S_values[idx], I_values[idx], R_values[idx])
end

# 显示求解器信息
println()
println("求解器信息：")
println("使用DP5方法（Dormand-Prince 5阶方法，即ode45）")
println("相对误差容限：1e-8")
println("绝对误差容限：1e-8")
println("总步数：$(length(sol.t)) (刚体运动)")
println("总步数：$(length(sol_sir.t)) (SIR模型)")

# 保存数据到文件以便在Syslab中绘图
println()
println("正在保存数据到文件...")

# 保存刚体运动数据（使用英文列名）
open("rigid_body_data.txt", "w") do io
    println(io, "time,y1,y2,y3")
    for i in 1:length(t_values)
        println(io, "$(t_values[i]),$(y1_values[i]),$(y2_values[i]),$(y3_values[i])")
    end
end

# 保存SIR模型数据（使用英文列名）
open("SIR_data.txt", "w") do io
    println(io, "time,susceptible,infected,recovered")
    for i in 1:length(t_sir_values)
        println(io, "$(t_sir_values[i]),$(S_values[i]),$(I_values[i]),$(R_values[i])")
    end
end

println("数据已保存到文件：")
println("- rigid_body_data.txt")
println("- SIR_data.txt")
println()
println("在Syslab中绘制图形的代码：")
println("// 读取数据")
println("data = readtable('rigid_body_data.txt');")
println("t = data.time;")
println("y1 = data.y1;")
println("y2 = data.y2;")
println("y3 = data.y3;")
println()
println("// 绘制图形")
println("figure;")
println("plot(t, y1, 'b-o', 'LineWidth', 2, 'MarkerSize', 6, 'DisplayName', 'y₁');")
println("hold on;")
println("plot(t, y2, 'r-o', 'LineWidth', 2, 'MarkerSize', 6, 'DisplayName', 'y₂');")
println("plot(t, y3, 'y-o', 'LineWidth', 2, 'MarkerSize', 6, 'DisplayName', 'y₃');")
println("title('刚性体运动方程的数值解', 'FontSize', 14, 'FontWeight', 'bold');")
println("xlabel('时间', 'FontSize', 12);")
println("ylabel('状态变量', 'FontSize', 12);")
println("grid on;")
println("legend('Location', 'southwest');")
println("xlim([0, 12]);")
println("ylim([-1.2, 1.2]);")