如何将MATLAB脚本转换为返回x与u(x,t)的函数？

作为刚接触MATLAB的新手，这种脚本转函数的坑我当年也踩过！首先我帮你补全了那段没写完的初始条件（看起来是波动方程常用的高斯脉冲测试用例），然后整理出正确的函数写法，再说说你大概率踩的那些导致绘图不对的坑。

先还原你原脚本的完整逻辑（补全截断部分） #

原脚本应该是用显式差分法求解一维波动方程，完整代码大概是这样：

xmin=0; xmax=1; N=100; dt=0.009; t=0; tmax=0.5; v=1;
% 离散化计算域
dx=(xmax-xmin)/N; 
x=xmin-dx : dx : xmax+dx; % 加虚拟边界点用于差分计算
% 初始条件：中心在0.5的高斯脉冲
u0=exp(-200*(x-0.5).^2); 
u_prev = u0; % 存储t-dt时刻的解
% 计算第一个时间步（波动方程需要前两个时间层的解）
r = v*dt/dx; % Courant数，稳定条件要求r<=1
u_current = u_prev + 0.5*r^2 * [0, u_prev(3:end)-2*u_prev(2:end-1)+u_prev(1:end-2), 0];
t = t + dt;
% 时间循环迭代到tmax
while t < tmax
    u_next = 2*u_current - u_prev + r^2 * [0, u_current(3:end)-2*u_current(2:end-1)+u_current(1:end-2), 0];
    % 更新时间层变量
    u_prev = u_current;
    u_current = u_next;
    t = t + dt;
end
% 绘图
plot(x, u_current);
xlabel('x'); ylabel('u(x,t)');
title(['t = ', num2str(t)]);

转换后的可返回x和u(x,t)的函数 #

我把这段脚本改成了标准的MATLAB函数，完全保留原脚本的计算逻辑，同时保证变量正确传递出来：

function [x, u_final] = solve_wave_equation(tmax)
    % 设置默认参数：如果调用时不传入tmax，就用0.5
    if nargin < 1
        tmax = 0.5;
    end
    
    % 原脚本里的固定参数（你也可以改成输入参数，让函数更灵活）
    xmin = 0;
    xmax = 1;
    N = 100;
    dt = 0.009;
    v = 1;
    t = 0;
    
    % 离散化空间域，保留虚拟边界点
    dx = (xmax - xmin)/N;
    x = xmin - dx : dx : xmax + dx;
    
    % 初始条件：和原脚本一致的高斯脉冲
    u0 = exp(-200*(x - 0.5).^2);
    u_prev = u0; % t-dt时刻的解
    
    % 计算第一个时间步
    r = v*dt/dx;
    u_current = u_prev + 0.5*r^2 * [0, u_prev(3:end)-2*u_prev(2:end-1)+u_prev(1:end-2), 0];
    t = t + dt;
    
    % 时间循环迭代到目标时刻
    while t < tmax
        u_next = 2*u_current - u_prev + r^2 * [0, u_current(3:end)-2*u_current(2:end-1)+u_current(1:end-2), 0];
        
        % 更新时间层
        u_prev = u_current;
        u_current = u_next;
        t = t + dt;
    end
    
    % 最终返回的时刻解
    u_final = u_current;
end

正确调用函数并绘图的方法 #

你只需要在命令行或者另一个脚本里写这段代码，就能得到和原脚本一模一样的绘图效果：

% 调用函数，获取x坐标和t=0.5时的u(x,t)
[x, u] = solve_wave_equation(0.5);
% 绘图
plot(x, u);
xlabel('x'); ylabel('u(x,t)');
title('Wave Equation Solution at t=0.5');
grid on;

你之前绘图效果不对的大概率原因 #

1. 变量作用域踩坑：脚本里的变量都是全局的，但函数里的变量默认是局部的——如果你的函数没有正确返回x和u_final，或者绘图时误用了函数外的未定义变量，肯定会出问题。
2. 边界点丢失：原脚本里的x包含了xmin-dx和xmax+dx这两个虚拟边界点，用来处理差分的边界条件。如果你的函数里不小心删掉了这两个点，差分计算就会出错，波形直接变形。
3. 时间循环逻辑错了：比如把循环条件写成t <= tmax，导致多跑了几个时间步，得到的不是你想要的时刻的解；或者循环根本没执行足够次数，解还停留在初始时刻。
4. 初始条件补错了：你原代码里的u0=exp(-200...没写完，如果补成了别的形式，自然绘图效果和原脚本不一样。

给新手的额外小提示 #

如果想让函数更灵活，比如可以自己设置x范围、时间步这些参数，你可以把函数改成支持多输入参数的版本：

function [x, u_final] = solve_wave_equation(xmin, xmax, N, dt, tmax, v)
    % 给每个参数设置默认值，调用时可以缺省后面的参数
    if nargin < 6, v = 1; end
    if nargin < 5, tmax = 0.5; end
    if nargin < 4, dt = 0.009; end
    if nargin < 3, N = 100; end
    if nargin < 2, xmax = 1; end
    if nargin < 1, xmin = 0; end
    
    % 后面的计算逻辑和之前的函数完全一致
    t = 0;
    dx = (xmax - xmin)/N;
    x = xmin - dx : dx : xmax + dx;
    u0 = exp(-200*(x - 0.5).^2);
    u_prev = u0;
    r = v*dt/dx;
    u_current = u_prev + 0.5*r^2 * [0, u_prev(3:end)-2*u_prev(2:end-1)+u_prev(1:end-2), 0];
    t = t + dt;
    
    while t < tmax
        u_next = 2*u_current - u_prev + r^2 * [0, u_current(3:end)-2*u_current(2:end-1)+u_current(1:end-2), 0];
        u_prev = u_current;
        u_current = u_next;
        t = t + dt;
    end
    
    u_final = u_current;
end

这样你可以根据需求灵活调整参数，比如solve_wave_equation(0,2,200,0.005,1,1.5)就能求解x范围0-2、200个网格点、时间步0.005、到t=1、波速1.5的波动方程。

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