checkjac.m

来自「电力系统分析计算程序」· M 代码 · 共 37 行

numjacs
fm_call('i')

disp(' ')

disp(['Fx abs err = ',num2str(max(max(abs(DAE.Fx-Fx))))])
disp(['Fy abs err = ',num2str(max(max(abs(DAE.Fy-Fy))))])
disp(['Gx abs err = ',num2str(max(max(abs(DAE.Gx-Gx))))])
disp(['Gy abs err = ',num2str(max(max(abs(DAE.Gy-Gy))))])

disp(' ')

if DAE.n 
  [i,j] = find(abs(DAE.Fx) <= Settings.dyntol);
  Fxtmp = DAE.Fx + sparse(i,j,1,DAE.n,DAE.n);
  
  [i,j] = find(abs(DAE.Fy) <= Settings.dyntol);
  Fytmp = DAE.Fy + sparse(i,j,1,DAE.n,DAE.m);
  
  [i,j] = find(abs(DAE.Gx) <= Settings.dyntol);
  Gxtmp = DAE.Gx + sparse(i,j,1,DAE.m,DAE.n);
else
  Fxtmp = 1;
  Fytmp = ones(1,DAE.m);
  Gxtmp = ones(DAE.m,1);
end

[i,j] = find(abs(DAE.Gy) <= Settings.dyntol);
Gytmp = DAE.Gy + sparse(i,j,1,DAE.m,DAE.m);

disp(['Fx rel err = ',num2str(max(max(abs((DAE.Fx-Fx)./Fxtmp))))])
disp(['Fy rel err = ',num2str(max(max(abs((DAE.Fy-Fy)./Fytmp))))])
disp(['Gx rel err = ',num2str(max(max(abs((DAE.Gx-Gx)./Gxtmp))))])
disp(['Gy rel err = ',num2str(max(max(abs((DAE.Gy-Gy)./Gytmp))))])

disp(' ')