I am looking for finding the velocity distribution using Principle of Maximum entropy when applied to Tsallis entropy.Tsallis entropy is defined as:$$S_{T} = \frac{1}{q-1}\left(1-\int{f(p)^qdp}\right)$$I have the following constraints:$$\int{f(p)dp}=1\\\int{pf(p)dp}=0\\\int{p^2f(p)dp}=1\\\int{p^3f(p)dp}=0\\$$Using these constraints and incorporating Lagrange multipliers, the equation of ${S_{T}}$ becomes:$$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)$$Taking variation of ${S_{T}}$ with respect to ${f(p)}$ and equating it to 0, the distribution function becomes:$$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}}$$The entire problem then boils down to solving for the Lagrange multiplier. However, I am unable to solve this problem using Newton-Raphson. Can anyone give me some pointers as to how to proceed with solving this problem. My MATLAB Newton Raphson code incorporating first three constraints is shown below. Any help would be appreciated:%Here we use Newton_Solverclear all;clc;%Initial Trialslambda(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 indexwhile(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
Constrained Minimization of Tsallis Entropy
algorithms