소개

이번에는 일본어 버전이 나오지 않았던 아래의 문서를 일본어로 번역해 나가겠습니다. (나의 작은 말이 더 많을지도 모른다. 별로 축어 번역하는 기분은 없습니다.)
htps : // jp. 마 t 후 rks. 이 m/헤일프/파랏ぇㅇㅇㅇㅇㅁㅁㅁㅁㅁㅁ HTML

사례: 로렌츠계

이번에는 이하의 상미분 방정식으로 표현되는 로렌츠계에서의 파라미터 스윕을 실행합니다. $\sigma,\rho,\beta$ 는 해당 시스템의 매개 변수를 나타냅니다.

\begin{align}
\frac{\mathrm{d}}{\mathrm{d}t}x &= \sigma(y-z) \\
\frac{\mathrm{d}}{\mathrm{d}t}y &= x(\sigma - z) -y \\
\frac{\mathrm{d}}{\mathrm{d}t}z &= xy - \beta x
\end{align}

절각이므로, Loretz(1963) 1 가 상정한 계를 실장해 보겠습니다. 솔버는 ode45를 사용했습니다. (평상시부터 ode45 로 해 보고, 좋다면 요시!같이 하고 있습니다.여기엔 공부 부족...) a 의 각 행이 위의 수식의 $x,y,z$ 에 대응 하고 있습니다.

Lorentz_System_1963

sigma = 10;
beta = 8/3;
rho = 24.74;
f = @(t,a) [sigma*(-a(1) + a(2)); a(1)*(rho - a(3)) - a(2); a(1)*a(2) - beta*a(3) ];
[t,a] = ode45(f,[0 100],[1 1 1]); % [x0 y0 z0] = [1 1 1]
plot3(a(:,1),a(:,2),a(:,3))
xlabel('x'); ylabel('y'); zlabel('z')
title('Lorentz(1963)')

원하는 그림을 얻었습니다.

매개 변수 스윕을 실행해보십시오.

이제 매개 변수 스윕에 대한 함수를 만들어 보겠습니다. 매개 변수의 인덱스 너비, 매개 변수 및 대기열을 입력 인수로 사용하고 최종 도달 지점의 z 좌표를 result에 저장합니다. 인덱스 폭을 지정하는 의미는 나중에.

parameterSweep.m

function results = parameterSweep(first,last,sigma,rho,beta,Q)
    results = zeros(last-first,1);
    for ii = first:last-1
        lorenzSystem = @(t,a) [sigma(ii)*(a(2) - a(1)); a(1)*(rho(ii) - a(3)) - a(2); a(1)*a(2) - beta*a(3)];
        [t,a] = ode45(lorenzSystem,[0 100],[1 1 1]);
        result = a(end,3); % z 座標を返す。
        send(Q,[ii,result]); % キュー (DataQueue Object) を用いて，ワーカーからクライアントに結果を送る。
        results(ii-first+1) = result;
    end
end

매개 변수 스윕을 실행하기 전에 좌표축도 설정하십시오. 이번에는 $\sigma,\rho$ 를 변화시켜 보고 어떤 변화가 일어나는지를 시각화해 보겠습니다.

createParameterGrid.m

gridSize = 40;
sigma = linspace(5, 45, gridSize);
rho = linspace(50, 100, gridSize);
beta = 8/3;
[rho,sigma] = meshgrid(rho,sigma);
figure('Visible',true);
surface = surf(rho,sigma,NaN(size(sigma))); % z 座標を NaN にしておく。
xlabel('\rho','Interpreter','Tex'); ylabel('\sigma','Interpreter','Tex')

좌표축의 설정을 할 수 있었으므로, 이 위에 플롯 해 갑니다. 먼저 parallel pool을 시작하고 DataQueue Object를 만듭니다.

parpool_and_queue.m

parpool; % parallel pool 起動
Q = parallel.pool.DataQueue
afterEach(Q,@(data) updatePlot(surface,data));

updatePlot 함수는 다음과 같습니다.

updatePlot.m

function updatePlot(surface,data)
    surface.ZData(data(1)) = data(2); % 上で NaN にしていた z 座標に代入する。
    drawnow('limitrate'); % 更新のレートを 20 frames /sec にする。
end

위의 createParameter 함수에서 인덱스 너비를 했는데, 이는 작업자에게로드 할 매개 변수를 분할하기 때문입니다. parfeval 작업자에 대한로드를 분할하는 것이 더 효율적으로 작동합니다. 분할의 기준에 관해서는, 문서에서는 「하나의 태스크가, 잡 스케줄링에 의한 오버헤드보다 충분히 큰 계산 시간을 소비하는 사이즈」로, 「태스크수가, 모든 워커에 할당할 수 있을 뿐」이라고 설명되어 있습니다. (엄청 의역하고 있습니다만, 분산 처리의 결정이군요.)
그런 다음 $\sigma,\rho$를 다음과 같이 분할하고 매개 변수 스윕을 수행합니다.

parameterSweepCommand.m

step = 100;
partitions = [1:step:numel(sigma), numel(sigma)+1];

f  = parallel.FevalFuture; % ここ仕様変更(R2019b 時点)があったみたいで，公式から変更していますが，ただの事前割り当てです。
for ii = 1:numel(partitions)-1
    f(ii) = parfeval(@parameterSweep,1,partitions(ii),partitions(ii+1),sigma,rho,beta,Q);
end

출력할 수 있었습니다. 작업자 수 8에서 실행 중입니다. 잘 보면, 4x2=8씩 플롯이 진행되고 있는 것을 알 수 있군요. (MATLAB에서 플롯을 gif로 출력하는 기술을 착용하지 않으므로 캡처 (쓴웃음))wait(f);를 입력하면 f 타이머가 실행을 완료 할 때까지 명령을 차단할 수 있으므로 위의 parfeval 결과에 의존하는 명령 오류를 방지 할 수 있습니다.

LorentzContour.m

results = reshape(fetchOutputs(f),gridSize,[]);
contourf(rho,sigma,results)
xlabel('\rho','Interpreter','Tex')
ylabel('\sigma','Interpreter','Tex')

글쎄, 이것이 무엇을 시사하는지 상쾌합니다. 「역시 카오스군요, 초기 조건 움직인 것은 아니지만」라고 하는 것으로 좋을까. 이 다큐멘테이션이 일본어 번역되지 않는 이유를 어쩐지 알 수 있다고 생각한다.

요약

parfeval 는 사전에 입력 인수를 분할해 주면 좋을 것이다.



E. N. Lorenz, “Deterministic Nonperiodic Flow,” Journal of Atmospheric Sciences, Vol. 20, No. 2, 1963, pp. 130-141.  ↩