当应用于 Tsallis 熵时，我正在寻找使用最大熵原理找到速度分布。

Tsallis 熵定义为：

$$S_{T} = \frac{1}{q-1}\left(1-\int{f(p)^qdp}\right)$$

我有以下约束：

$$\int{f(p)dp}=1\\ \int{pf(p)dp}=0\\ \int{p^2f(p)dp}=1\\ \int{p^3f(p)dp}=0\\ $$

使用这些约束并结合拉格朗日乘子，的方程变为：${S_{T}}$

$$S_{T} = \frac{1}{q-1}\left(1-\int{f(p)^qdp}\right) + \lambda_0 \left(\int{f(p)dp}-1\right) + \lambda_1 \left(\int{pf(p)dp}\right) + \lambda_2 \left(\int{p^2f(p)dp-1}\right)+ \lambda_3 \left(\int{p^3f(p)dp}\right)$$

取相对于的变化并将其等同于 0，分布函数变为： ${S_{T}}$${f(p)}$

$$f(p)=\left[\left(\frac{q-1}{q}\right)\left(\lambda_0+\lambda_1p+\lambda_2p^2+\lambda_3p^3\right)\right]^{\frac{1}{q-1}}$$

然后整个问题归结为求解拉格朗日乘数。但是，我无法使用 Newton-Raphson 解决这个问题。谁能给我一些关于如何继续解决这个问题的指示。我的包含前三个约束的 MATLAB Newton Raphson 代码如下所示。任何帮助，将不胜感激：

%Here we use Newton_Solver

clear all;
clc;


%Initial Trials
lambda(1,1) = -1.5;
lambda(2,1) = 0.0;
lambda(3,1) = 0.09;
lambda(4,1) = 0.09;

ll = -5;
ul = 5;
steps = 10000;
iter = 1500;
i = 0;
q = 1.5; %Tsallis q index
while(i<iter)
    func_val = zeros(4,1);
    Jacobian_mat = zeros(4,4);

    avg_v0 = 0.0;
    avg_v1 = 0.0;
    avg_v2 = 0.0;
    avg_v3 = 0.0;
    avg_v4 = 0.0;
    avg_v5 = 0.0;
    avg_v6 = 0.0;    
    for(j=1:steps)
        p = ll+(ul-ll)*j/steps;
        temp = (lambda(1) + lambda(2)*p + lambda(3)*p*p + lambda(4)*p*p*p)*(q-1)/q;
        temp = temp^(1/(q-1));
        temp = temp*(ul-ll)/steps;

        func_val(1) = func_val(1) + temp;
        func_val(2) = func_val(2) + temp*p;
        func_val(3) = func_val(3) + temp*p*p;
        func_val(4) = func_val(4) + temp*p*p*p;

        temp1 = lambda(1) + lambda(2)*p + lambda(3)*p*p + lambda(4)*p*p*p;
        temp1 = temp1^((2-q)/(q-1));
        temp1 = temp1*((q-1)/q)^(1/(q-1))/(q-1);
        temp1 = temp1*(ul-ll)/steps;

        avg_v0 = avg_v0 + temp1;
        avg_v1 = avg_v1 + temp1*p;
        avg_v2 = avg_v2 + temp1*p*p;
        avg_v3 = avg_v3 + temp1*p*p*p;
        avg_v4 = avg_v4 + temp1*p*p*p*p;
        avg_v5 = avg_v5 + temp1*p*p*p*p*p;
        avg_v6 = avg_v6 + temp1*p*p*p*p*p*p;

    end

    Jacobian_mat(1,1) = avg_v0;    
    Jacobian_mat(1,2) = avg_v1;   
    Jacobian_mat(1,3) = avg_v2;
    Jacobian_mat(1,4) = avg_v3;

    Jacobian_mat(2,1) = avg_v1;    
    Jacobian_mat(2,2) = avg_v2;
    Jacobian_mat(2,3) = avg_v3;
    Jacobian_mat(2,4) = avg_v5;

    Jacobian_mat(3,1) = avg_v2;    
    Jacobian_mat(3,2) = avg_v3;
    Jacobian_mat(3,3) = avg_v4;
    Jacobian_mat(3,4) = avg_v5;

    Jacobian_mat(4,1) = avg_v3;    
    Jacobian_mat(4,2) = avg_v4;
    Jacobian_mat(4,3) = avg_v5;
    Jacobian_mat(4,4) = avg_v6;

    func_val(1) = func_val(1) - 1.0;
    func_val(2) = func_val(2) - 0.0;
    func_val(3) = func_val(3) - 1.0;
    func_val(4) = func_val(4) - 0.0;

    Jacobian_inv = Jacobian_mat^(-1);
    lambda = lambda - Jacobian_inv*func_val*1.0;

    Disp_val = ['Step is: ',num2str(i), ' Lambda 1 values:', num2str(lambda(1,1))];
    Disp_val1 = ['Step is: ',num2str(i), ' Lambda 2 values:', num2str(lambda(2,1))];
    Disp_val2 = ['Step is: ',num2str(i), ' Lambda 3 values:', num2str(lambda(3,1))];
    Disp_val3 = ['Step is: ',num2str(i), ' Lambda 4 values:', num2str(lambda(4,1))];
    Disp_val4 = ['Step is: ',num2str(i), ' AVG 1 values:', num2str(func_val(1) + 1.0)];
    Disp_val5 = ['Step is: ',num2str(i), ' AVG 2 values:', num2str(func_val(2) + 0.0)];
    Disp_val6 = ['Step is: ',num2str(i), ' AVG 3 values:', num2str(func_val(3) + 1.0)];
    Disp_val7 = ['Step is: ',num2str(i), ' AVG 4 values:', num2str(func_val(4) + 0.0)];

    disp(Disp_val);
    disp(Disp_val1);
    disp(Disp_val2);
    disp(Disp_val3);
    disp(Disp_val4);
    disp(Disp_val5);
    disp(Disp_val6);
    disp(Disp_val7);
    disp(' ');
    i = i+1;
end