Gamma 函数的特殊值计算

gamma 函数在部分点的实现有些问题。例如 gamma(-4) Matlab 选择的是 +Inf，目前北太会给出错误。可以和 Matlab 行为保持一致，或者直接返回 NaN。这样对画图方便一些。目前的实现要画图，就只能手动分段。计算 1/gamma(x) 也变成不连续的了。即 Matlab 可以执行以下代码并画图

gamma([-5 -4 -3 -2 -1 0 5])

x = -5:0.01:5;
plot(x, gamma(x))

 Matlab R2023b 输出

>> gamma([-5 -4 -3 -2 -1 0 5])

ans = 

    Inf   Inf   Inf   Inf   Inf   Inf    24

北太 3.1.0 目前会报定义域错误

>> x = -5:0.01:5;
>> plot(x, gamma(x));
错误使用函数 gamma
domain error
程序执行中显示有错误信息，请反馈给开发团队。

目前要绘制  gamma 函数只能手动分段绘图：代码如下：

% 生成开区间 (start:stop)
function range=openRange(start, step, stop)
    range = (start+eps(start)):step:(stop-eps(stop));
end
% 绘制 gamma 函数
function gamma_plot()
    step = 0.01;
    x1 = [
        openRange(-5.0, step, -4.0)
        openRange(-4.0, step, -3.0)
        openRange(-3.0, step, -2.0)
        openRange(-2.0, step, -1.0)
        openRange(-1.0, step, -0.0)
    ]';
    x2 = [ openRange(0.0, step, 5.0) ]';
    
    % draw
    plot(...
        x1,gamma(x1),...
        x2,gamma(x2),...
        'LineWidth',2,...
    );
    grid on;
    xlim([-5, 5]);
    ylim([-10, 10]);
    title('Gamma(x) Line Plot')
    xlabel('x');
    ylabel('Gamma(x)');
end

对于高数里的极限、导数的计算有吗？积分有的

想尝试将高数里的使用MATLAB计算函数极限、导数、积分和求解微分方程转换为使用baltamatica，对于学习高数的学生而言，基本的计算能够解决就可以将此软件推广起来优点一：支持中文变量、中文符号-这个就很符合日常的习惯优点二：软件相比MATLAB可是要小太多啦

abs命令或者变量的类型发生变化导致,无法求值

写了一个自适应 Simpson 求积公式的代码, 运行过程中出现了 abs命令或者变量的类型发生变化导致,无法求值 的错误, 如图所示我写的 adapsimp 的代码如下

function [s, err] = adapsimp(func, a, b, tol)
    s = comsimp(func, a, b, 2);
    c = (a + b) / 2;
    s1 = comsimp(func, a, c, 2);
    s2 = comsimp(func, c, b, 2);
    s12 = s1 + s2;
    err = abs(s12 - s) / 15;
    if err < tol
        s = s12;
    else
        [s1, err1] = adapsimp(func, a, c, tol/2);
        [s2, err2] = adapsimp(func, c, b, tol/2);
        s = s1 + s2;
        err = err1 + err2;
    end
end

里面用到了 comsimp 函数, 是这样写的

function s = comsimp(func, a, b, n)
    h = (b - a) / n;
    s0 = func(a) + func(b);
    s1 = 0;     % summation of f(x_{2k-1})
    s2 = 0;     % summation of f(x_{2k})
    for k = 1:n-1
        x = a + k * h;
        if rem(k , 2) == 0
            s2 = s2 + func(x);
        else
            s1 = s1 + func(x);
        end
    end
    s = h * (s0 + 4 * s1 + 2 * s2) / 3;
end

这部分代码在 octave 上运行是没有问题的用的版本是2.2.0最新版的.

运行这段在MATLAB中写的迎风法解双曲型PDE时稳定闪退

h=0.01;
r=0.5;
tao=r*h;
x=0:h:5;
t=0:tao:10;
length_t=length(t);
length_x=length(x);
U=zeros(length_t,length_x);
length_in=length_x-1;
A=diag((1-r)*ones(length_in,1))+r*diag(ones(length_in-1,1),-1);
[X,T]=meshgrid(x,t);
F=2*T.*sin(X)+T.^2.*cos(X);
for i=2:length(t)
    U(i,2:end)=A*U(i-1,2:end)'+tao*F(i,2:end)';
end
U_acc=T.^2.*sin(X);
figure(1)
subplot(1,2,1)
mesh(x,t,U)
title("解曲面")
subplot(1,2,2)
mesh(x,t,U_acc-U)
title("误差曲面")

发现是FOR循环内的矩阵转置处出现了问题，但是无法清晰定位具体发生了什么导致的问题，希望开发组能解决。