求助：将求解递推多项式的Fortran程序转换为Matlab程序

我编写了如下Fortran程序，用于计算满足特定递推关系的多项式系数：

implicit real*16 (a-h,o-z)
    implicit integer*8 (i-n)
    
    parameter (nnmax=130000, nbarmax = nnmax/2,ipmax =5)
    parameter (nulatmax=5000)
    dimension doef(0:nbarmax)
    dimension coef(0:ipmax,0:nnmax)
    dimension aux(-1:nnmax)
    
    ip = 2
    anzn = 3
    
    coef = 0.q0 !-->zera todos os coeficientes
    
    alamb = 1.q0
    
    !Initial conditions
    do i1 =0,ip
       coef(i1,0) = 1.q0
    end do
!   write(*,*)'linha 271'
    aux(-1) = 0.q0
!   iteracao
    
    M = 30
    do nn =1, M
      nnbar = (nn+ip)/(ip+1)
      natual =mod(nn +2*(ip+1)-1,ip+1)
      nanter =mod(nn +2*(ip+1)-2,ip+1)
      nfirst =mod(nn +2*(ip+1)-(ip+2),ip+1)
      do i2=0,nnbar
        aux(i2) = coef(natual,i2)
      end do
      
      do i2=0,nnbar
        coef(natual,i2) = coef(nanter,i2) -(alamb**anzn)*aux(i2-1)
        doef(i2) = coef(natual,i2)
        write(511,*) nn,doef(i2)
      enddo
    enddo
          
endprogram

该程序对应的递推关系为：
$$P_{m+3}(x) = P_{m+2}(x) - x \cdot P_m(x) \quad (m \geq 0)$$
初始条件：

• $P_0(x) = 1$
• $P_1(x) = 1 - x$
• $P_2(x) = 1 - 2x$

我需要将上述计算多项式系数的Fortran例程改写为Matlab程序，并用Matlab的roots(X)函数求解这些多项式的根（其中X为多项式系数列向量）。

Matlab程序实现 #

以下是对应功能的Matlab代码，包含多项式系数计算和根求解两部分：

% 设置参数
ip = 2;
anzn = 3;
alamb = 1.0;
M = 30; % 递推次数，对应计算到P_{M+2}(x)

% 初始化系数矩阵：coef(i,j)代表P_{i-1}(x)中x^{j-1}项的系数（适配Matlab1-based索引）
max_degree = floor(M / (ip+1)); % 多项式最高次数
coef = zeros(ip+1, max_degree + 1); 

% 设置初始条件
coef(1, 1) = 1; % P_0(x)=1，x^0系数为1
coef(2, 1) = 1; coef(2, 2) = -1; % P_1(x)=1-x
coef(3, 1) = 1; coef(3, 2) = -2; % P_2(x)=1-2x

aux = zeros(1, max_degree + 2); % 辅助数组，预留i2-1的索引空间
aux(1) = 0; % 对应Fortran的aux(-1)=0，Matlab索引偏移适配

% 递推计算多项式系数
for nn = 1:M
    nnbar = floor((nn + ip) / (ip + 1));
    natual = mod(nn + 2*(ip+1) - 1, ip+1);
    natual = natual == 0 ? ip+1 : natual; % 转换为1-based索引
    nanter = mod(nn + 2*(ip+1) - 2, ip+1);
    nanter = nanter == 0 ? ip+1 : nanter;
    
    % 保存当前natual行的系数到aux
    for i2 = 0:nnbar
        aux(i2 + 2) = coef(natual, i2 + 1);
    end
    
    % 更新coef矩阵
    for i2 = 0:nnbar
        aux_prev = aux(i2 + 1); % 对应Fortran的aux(i2-1)
        coef(natual, i2 + 1) = coef(nanter, i2 + 1) - (alamb^anzn)*aux_prev;
    end
end

% 提取各多项式并求解根
for k = 0:ip
    % 获取P_k(x)的系数（x^0到最高次）
    poly_coeff = coef(k+1, :);
    % 去除末尾无效零系数
    poly_coeff = poly_coeff(1:find(poly_coeff~=0, 1, 'last'));
    % 反转顺序适配roots函数（需要最高次到常数项的顺序）
    poly_coeff_rev = flip(poly_coeff);
    
    fprintf('===== 多项式P_%d(x) =====\n', k);
    fprintf('系数（x^0到x^n）：');
    disp(poly_coeff);
    if length(poly_coeff) > 1
        roots_p = roots(poly_coeff_rev);
        fprintf('根：\n');
        disp(roots_p);
    else
        fprintf('该多项式为常数，无实根或复根\n');
    end
    fprintf('\n');
end

代码说明 #

• 索引适配：Matlab数组采用1-based索引，因此对原Fortran的0/负索引做了偏移处理，保证递推逻辑完全一致。
• 初始条件映射：严格按照题目给定的$P_0, P_1, P_2$初始化系数矩阵，确保计算起点正确。
• 递推逻辑复刻：完全复现Fortran中的循环和系数更新规则，保证计算结果与原程序一致。
• 根求解适配：roots函数要求系数按最高次到常数项排列，因此对提取的系数做反转处理，并去除末尾无效零值。

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