如何运行给定的微分方程MATLAB程序并替换自定义方程

1 内置微分方程位置说明 #

你提供的代码采用有限差分法求解二阶线性非齐次常微分方程，内置方程的定义分为三个部分：

• 方程系数定义：开头a=1;b=-2;c=2;对应通用二阶方程格式 a*y'' + b*y' + c*y = g(t) 里的系数，当前内置方程的系数对应 y'' - 2y' + 2y = g(t)
• 非齐次项定义：循环内t=k*dt;之后的g=exp(-t);就是非齐次项g(t)的表达式，当前内置为指数衰减函数
• 差分递推核心：循环内的y=((2+b*dt-c*dt^2)*y0+g*dt^2-y1)/(1+b*dt);是微分方程离散后的递推计算式，是整个求解逻辑的核心

2 自定义微分方程替换方法 #

根据你要求解的方程类型，按以下步骤修改即可：

• 如果你要解的是二阶线性非齐次常微分方程，直接替换开头的a、b、c参数为你的方程对应系数，再替换循环内g=exp(-t);为你自定义的非齐次项表达式即可
• 如果你要解的是其他类型微分方程，替换上述差分递推核心行即可，注意变量定义：y0为上一时刻的y值，y1为上上个时刻的y值，dt为时间步长，计算得到的y为当前时刻的y值
• 可调通用参数：修改dt=.001;可以调整时间步长精度，修改最外层for k=1:4000;的4000可以调整总模拟时长，修改pause(.1)可以调整动画的刷新速度

3 优化后带注释的完整代码 #

% 清空工作区与窗口
clc 
close all 
clear all

%% 微分方程参数定义 对应通用格式 a*y'' + b*y' + c*y = g(t)
a = 1;
b = -2;
c = 2;
% 系数归一化
b = b/a;
c = c/a;

%% 仿真与绘图参数定义
dt = 0.001;       % 时间步长
tap = 1;          % 纵轴范围基准值
btsa = tap;       % 纵轴上限
btsb = -tap;      % 纵轴下限
n = 40;           % 单摆绳锯齿采样点数
total_step = 4000;% 总仿真步数

%% 初始条件定义
y0 = 0;           % t=0时刻的y值
y1 = 1;           % t=-dt时刻的y值，用于递推启动
y1 = y0 - y1*dt;

%% 预生成单摆绳的x坐标（固定锯齿形状）
for k = 1:n
    px(k) = (-1)^k * 0.1;
end
px = [0 0 px 0 0]; % 补全两端点

g0 = 0; % 上一时刻的非齐次项值，用于绘制激励曲线
g = 0;  % 当前时刻的非齐次项值

%% 仿真主循环
for k = 1:total_step
    t = k*dt;
    
    % 非齐次项g(t)计算，可替换为自定义表达式
    g = exp(-t);
    g = g/a;
    
    % 核心差分递推计算当前时刻y值，自定义方程替换此处即可
    y = ((2 + b*dt - c*dt^2)*y0 + g*dt^2 - y1)/(1 + b*dt);
    
    %% 动画子图（2,2,1）绘制：单摆动画
    by = y + 0.3;
    py = linspace(by, tap-0.3, n+1) - (tap-0.3 - by)/(2*n);
    py(1) = py(1) + (tap-0.3 - by)/(2*n);
    py = [y py tap-0.3 tap];
    
    subplot(2,2,1)
    hold off
    plot([-2 2],[tap tap],'Color','k','LineWidth',4) % 悬挂固定梁
    hold on
    plot([0 4],[y y],'Color',[.8 .8 .8]); % 平衡位置参考线
    plot(px,py,'-r') % 单摆绳
    plot(0,y,'--o','MarkerFaceColor',[1 0 0],'MarkerEdgeColor',[0 0 0],'MarkerSize',20) % 摆球
    axis([-3 3 -tap+.1 tap+.1])
    
    %% 子图（2,2,2）绘制：求解得到的y(t)曲线
    subplot(2,2,2)
    plot([0 k*dt+10],[0 0],'-k') % 零点参考线
    hold on
    gam = plot([0 (k+1)*dt],[y y],'Color',[.8 .8 .8]); % 当前值参考线
    plot([k*dt (k+1)*dt],[y0 y],'-b') % 最新一段y(t)曲线
    
    % 自适应调整纵轴范围
    btsa = max(max(btsa,g),y);
    btsb = min(min(btsb,g),y);
    axis([0 max(10,k*dt)+.5 btsb-.1 btsa+.1])
    
    % 更新历史y值，为下一次递推做准备
    y1 = y0;
    y0 = y;
    
    %% 子图（2,2,4）绘制：激励项g(t)曲线
    subplot(2,2,4)
    plot([0 k*dt+10],[0 0],'-k') % 零点参考线
    hold on
    plot([k*dt (k+1)*dt],[g0 g],'-r') % 最新一段g(t)曲线
    axis([0 max(10,k*dt)+.5 btsb-.1 btsa+.1])
    
    % 更新历史g值
    g0 = g;
    
    % 控制动画刷新速度
    pause(0.1)
    delete(gam)
end

内容的提问来源于stack exchange，提问作者imaginary_eye