昆仑山:眼中无形心中有穴之穴人合一

夫君子之行,静以修身,俭以养德;非澹泊无以明志,非宁静无以致远。夫学须静也,才须学也;非学无以广才,非志无以成学。怠慢则不能励精,险躁则不能冶性。年与时驰,意与岁去,遂成枯落,多不接世。悲守穷庐,将复何及!

 

交互式动画演示经典的四阶龙格 - 库塔方法(RK4)求解常微分方程 (ODE) 的初值问题

提示词

我将使用 Python 开发一个交互式动画演示程序,通过经典的四阶龙格 - 库塔方法(RK4)求解常微分方程 (ODE) 的初值问题。这个程序将集成以下功能:

使用 matplotlib 实现动态可视化,展示 ODE 解的演化过程
通过 tkinter 构建直观的参数控制面板,允许用户实时调整:
初始条件参数
时间范围设置
RK4 算法的步长
ODE 方程的系数
程序布局设计为上下结构:
上半部分为动画演示区,动态展示 ODE 解曲线
下半部分为参数控制面板,包含滑动条、输入框和控制按钮
实时交互功能:
可以暂停 / 继续动画
调整参数后即时重新计算并更新显示
保存当前动画为视频文件

源代码


import tkinter as tk
from tkinter import ttk
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
from scipy.integrate import solve_ivp
import matplotlib
import warnings

# 设置中文字体和负号显示
try:
    matplotlib.rcParams['font.family'] = 'Microsoft YaHei'
except Exception:
    warnings.warn("系统未安装 'Microsoft YaHei',已使用默认字体。")
matplotlib.rcParams['axes.unicode_minus'] = False


class RK4App:
    """
    RK4方法求解常微分方程的图形化演示应用
    
    该应用使用Tkinter创建GUI界面,通过Matplotlib可视化展示四阶龙格-库塔法(RK4)
    求解常微分方程的过程。包含数值解和精确解对比、k值计算步骤可视化等功能。
    
    属性:
        root: Tkinter根窗口对象
        k1_val, k2_val, k3_val, k4_val: 用于显示当前k值的字符串变量
        fig: Matplotlib图形对象
        ax_num, ax_exact: 数值解和精确解子图
        canvas: 嵌入Tkinter的Matplotlib画布
        func_str: 微分方程字符串
        x0, y0: 初始条件
        x_end: 求解区间终点
        h: 步长
        steps: 总步数
        x_vals, y_vals: 数值解坐标点
        k_vals: 存储k值
        current_step: 当前计算步数
        sub_step: RK4子步骤(0-4)
    """
    def __init__(self, root):
        """
        初始化RK4演示应用
        
        参数:
            root: Tkinter根窗口对象
        """
        self.root = root
        self.root.title("RK4 求解常微分方程动画演示")
        # k 值显示变量
        self.k1_val = tk.StringVar(value="k1 = ")
        self.k2_val = tk.StringVar(value="k2 = ")
        self.k3_val = tk.StringVar(value="k3 = ")
        self.k4_val = tk.StringVar(value="k4 = ")

        # 创建图形区域
        self.fig, (self.ax_num, self.ax_exact) = plt.subplots(1, 2, figsize=(12, 5))

        self.canvas = FigureCanvasTkAgg(self.fig, master=root)
        self.canvas.get_tk_widget().pack(side=tk.TOP, fill=tk.BOTH, expand=True)

        # 参数
        self.func_str = tk.StringVar(value="y - x**2 + 1")
        self.x0 = tk.DoubleVar(value=0.0)
        self.y0 = tk.DoubleVar(value=0.5)
        self.x_end = tk.DoubleVar(value=4.0)
        self.h = tk.DoubleVar(value=0.2)
        self.steps = tk.IntVar(value=20)

        self.create_controls()
        self.reset_state()

    def create_controls(self):
        """
        创建GUI控制元素
        
        包括参数输入框、控制按钮、k值显示面板等
        """
        frame = tk.Frame(self.root)
        frame.pack(side=tk.BOTTOM, fill=tk.X)

        def labeled_entry(parent, label_text, var, width=8):
            """创建带标签的输入框"""
            tk.Label(parent, text=label_text, width=8).pack(side=tk.LEFT)
            tk.Entry(parent, textvariable=var, width=width).pack(side=tk.LEFT)

        # 第一行:微分方程输入
        row1 = tk.Frame(frame); row1.pack(fill=tk.X, padx=10, pady=2)
        tk.Label(row1, text="dy/dx = ", width=8).pack(side=tk.LEFT)
        tk.Entry(row1, textvariable=self.func_str, width=40).pack(side=tk.LEFT, expand=True, fill=tk.X)

        # 第二行:初始条件输入
        row2 = tk.Frame(frame); row2.pack(fill=tk.X, padx=10, pady=2)
        labeled_entry(row2, "x0", self.x0)
        labeled_entry(row2, "y0", self.y0)

        # 第三行:步长和步数输入
        row3 = tk.Frame(frame); row3.pack(fill=tk.X, padx=10, pady=2)
        labeled_entry(row3, "步长 h", self.h)
        labeled_entry(row3, "步数", self.steps)

        # 第四行:终点输入
        row4 = tk.Frame(frame); row4.pack(fill=tk.X, padx=10, pady=2)
        labeled_entry(row4, "x_end", self.x_end)

        # 按钮区域
        btn_frame = tk.Frame(frame); btn_frame.pack(pady=5)

        # k值显示面板
        k_frame = tk.LabelFrame(frame, text="当前步的 k 值", padx=10, pady=5)
        k_frame.pack(fill=tk.X, padx=10, pady=5)

        tk.Label(k_frame, textvariable=self.k1_val, fg='red').pack(anchor='w')
        tk.Label(k_frame, textvariable=self.k2_val, fg='green').pack(anchor='w')
        tk.Label(k_frame, textvariable=self.k3_val, fg='purple').pack(anchor='w')
        tk.Label(k_frame, textvariable=self.k4_val, fg='blue').pack(anchor='w')

        # 控制按钮
        self.start_btn = ttk.Button(btn_frame, text="开始演示", command=self.start)
        self.next_btn = ttk.Button(btn_frame, text="下一步", command=self.next_step, state=tk.DISABLED)
        self.reset_btn = ttk.Button(btn_frame, text="重置", command=self.reset, state=tk.DISABLED)
        for btn in [self.start_btn, self.next_btn, self.reset_btn]:
            btn.pack(side=tk.LEFT, padx=5)

        # 信息显示标签
        self.info_label = tk.Label(frame, text="请设置参数,点击开始演示", anchor='w')
        self.info_label.pack(fill=tk.X, padx=10)

    def reset_state(self):
        """重置计算状态和可视化数据"""
        self.x_vals = []
        self.y_vals = []
        self.k_vals = []
        self.k_arrows = []
        self.current_step = 0
        self.sub_step = 0
        self.global_error = 0.0

    def start(self):
        """开始演示:初始化计算并绘制初始图形"""
        # 解析微分方程函数
        try:
            expr = self.func_str.get()
            self.func = eval("lambda x,y: " + expr, {"__builtins__": None},
                             {"np": np, "sin": np.sin, "cos": np.cos, "exp": np.exp})
        except Exception as e:
            self.info_label.config(text=f"函数解析错误: {e}")
            return

        # 重置状态变量
        self.reset_state()
        self.x_vals = [self.x0.get()]
        self.y_vals = [self.y0.get()]
        self.h_val = self.h.get()
        self.total_steps = self.steps.get()

        # 计算精确解
        try:
            sol = solve_ivp(self.func, (self.x0.get(), self.x_end.get()), [self.y0.get()],
                            method='RK45', rtol=1e-8, atol=1e-8)
            self.x_exact, self.y_exact = sol.t, sol.y[0]
        except Exception as e:
            self.info_label.config(text=f"精确解计算失败: {e}")
            return

        # 初始化数值解图形
        self.ax_num.clear()
        self.ax_num.set_title("RK4 数值解")
        self.ax_num.set_xlabel("x")
        self.ax_num.set_ylabel("y")
        self.ax_num.grid(True)

        # 初始化精确解图形
        self.ax_exact.clear()
        self.ax_exact.set_title("精确解")
        self.ax_exact.set_xlabel("x")
        self.ax_exact.set_ylabel("y")
        self.ax_exact.grid(True)
        self.ax_exact.plot(self.x_exact, self.y_exact, 'g--', label="精确解")
        self.ax_exact.legend()

        # 绘制初始数值解点
        self.line, = self.ax_num.plot(self.x_vals, self.y_vals, 'b-', label="RK4 解")
        self.point, = self.ax_num.plot(self.x_vals, self.y_vals, 'ro', label="当前点")
        self.ax_num.legend()

        # 更新UI状态
        self.start_btn.config(state=tk.DISABLED)
        self.next_btn.config(state=tk.NORMAL)
        self.reset_btn.config(state=tk.NORMAL)
        self.info_label.config(text="开始演示 RK4,点击下一步查看每步计算过程。")
        self.canvas.draw()

    def next_step(self):
        """执行下一步计算并更新可视化"""
        # 检查是否已完成所有步骤
        if self.current_step >= self.total_steps:
            self.info_label.config(text="已完成全部步骤。")
            self.next_btn.config(state=tk.DISABLED)
            return

        # 准备当前步骤的RK4计算
        if self.sub_step == 0:
            self.prepare_rk4_step()

        # 清除上一步的箭头标注
        self.clear_arrows()
        
        # 显示当前子步骤
        self.show_substep()
        self.canvas.draw()

        # 更新子步骤计数器
        self.sub_step += 1
        if self.sub_step > 4:
            self.sub_step = 0

    def prepare_rk4_step(self):
        """计算当前步骤的四个k值"""
        self.x_n = self.x_vals[-1]
        self.y_n = self.y_vals[-1]
        h = self.h_val
        try:
            # 计算四个k值
            self.k1 = h * self.func(self.x_n, self.y_n)
            self.k2 = h * self.func(self.x_n + h/2, self.y_n + self.k1/2)
            self.k3 = h * self.func(self.x_n + h/2, self.y_n + self.k2/2)
            self.k4 = h * self.func(self.x_n + h, self.y_n + self.k3)
            
            # 计算下一步的y值
            self.y_next = self.y_n + (self.k1 + 2*self.k2 + 2*self.k3 + self.k4)/6
            self.x_next = self.x_n + h
        except Exception as e:
            self.info_label.config(text=f"RK4计算失败: {e}")
            self.next_btn.config(state=tk.DISABLED)

    def show_substep(self):
        """可视化当前子步骤的计算过程"""
        step = self.current_step + 1
        h = self.h_val

        # 根据子步骤显示不同的计算信息
        if self.sub_step == 0:
            self.add_arrow(self.x_n, self.y_n, self.k1, 'r', 'k1')
            info = f"步骤 {step}: 计算 k1\nk1 = {self.k1:.4f}"

        elif self.sub_step == 1:
            x_mid = self.x_n + h/2
            y_mid = self.y_n + self.k1/2
            self.add_arrow(x_mid, y_mid, self.k2, 'g', 'k2')
            info = f"步骤 {step}: 计算 k2\nk2 = {self.k2:.4f}"

        elif self.sub_step == 2:
            x_mid = self.x_n + h/2
            y_mid = self.y_n + self.k2/2
            self.add_arrow(x_mid, y_mid, self.k3, 'm', 'k3')
            info = f"步骤 {step}: 计算 k3\nk3 = {self.k3:.4f}"

        elif self.sub_step == 3:
            x_end = self.x_n + h
            y_end = self.y_n + self.k3
            self.add_arrow(x_end, y_end, self.k4, 'c', 'k4')
            info = f"步骤 {step}: 计算 k4\nk4 = {self.k4:.4f}"

        elif self.sub_step == 4:
            # 完成当前步骤,更新数值解
            self.x_vals.append(self.x_next)
            self.y_vals.append(self.y_next)
            self.line.set_data(self.x_vals, self.y_vals)
            self.point.set_data([self.x_next], [self.y_next])
            self.current_step += 1

            # 计算并更新误差
            idx = np.searchsorted(self.x_exact, self.x_next)
            idx = min(idx, len(self.x_exact)-1)
            y_exact = self.y_exact[idx]
            err = abs(self.y_next - y_exact)
            self.global_error = max(self.global_error, err)
            info = f"步骤 {step}: y_{step} = {self.y_next:.4f}, 误差 = {err:.2e}"

        # 更新UI显示
        self.info_label.config(text=info)
        self.k1_val.set(f"k1 = {self.k1:.5f}")
        self.k2_val.set(f"k2 = {self.k2:.5f}")
        self.k3_val.set(f"k3 = {self.k3:.5f}")
        self.k4_val.set(f"k4 = {self.k4:.5f}")

    def add_arrow(self, x, y, k, color, label):
        """
        在图形上添加表示导数的箭头
        
        参数:
            x, y: 箭头起点坐标
            k: 斜率值
            color: 箭头颜色
            label: 箭头标签
        """
        dx = self.h_val * 0.1
        dy = k * 0.1
        arrow = self.ax_num.arrow(x, y, dx, dy, head_width=0.05, head_length=0.05,
                                  fc=color, ec=color)
        self.ax_num.text(x + dx, y + dy, label, color=color)
        self.k_arrows.append(arrow)

    def clear_arrows(self):
        """清除所有k值箭头标注"""
        for arrow in self.k_arrows:
            arrow.remove()
        self.k_arrows.clear()

    def reset(self):
        """重置应用状态到初始界面"""
        # 重置内部状态
        self.reset_state()
        
        # 清除图形
        self.ax_num.clear()
        self.ax_num.set_title("RK4 数值解")
        self.ax_num.set_xlabel("x")
        self.ax_num.set_ylabel("y")
        self.ax_num.grid(True)

        self.ax_exact.clear()
        self.ax_exact.set_title("精确解")
        self.ax_exact.set_xlabel("x")
        self.ax_exact.set_ylabel("y")
        self.ax_exact.grid(True)

        # 重绘画布
        self.canvas.draw()
        
        # 更新UI状态
        self.start_btn.config(state=tk.NORMAL)
        self.next_btn.config(state=tk.DISABLED)
        self.reset_btn.config(state=tk.DISABLED)
        self.info_label.config(text="请设置参数,点击开始演示")

        # 重置k值显示
        self.k1_val.set("k1 = ")
        self.k2_val.set("k2 = ")
        self.k3_val.set("k3 = ")
        self.k4_val.set("k4 = ")

if __name__ == "__main__":
    root = tk.Tk()
    app = RK4App(root)
    root.mainloop()


posted on 2025-07-10 00:17  Indian_Mysore  阅读(19)  评论(0)    收藏  举报

导航