Programmering (forts), ett simuleringsexempel

Contents

Ett simuleringsexempel

% Fixa en 'startplåt'
T0 = zeros(5);
T0(2:4,1) = 80; T0(2:4,5) = 100;
T0(1,2:4) = 40; T0(5,2:4) = 20;
T0

T0 =

     0    40    40    40     0
    80     0     0     0   100
    80     0     0     0   100
    80     0     0     0   100
     0    20    20    20     0

% Beräkna temp i nästa tidssteg

[m,n] = size(T0);
T1 = T0;
for i = 2:m-1
    for j = 2:n-1
        T1(i,j) = (T0(i-1,j) + T0(i,j+1) + T0(i+1,j) + T0(i,j-1))/4;
    end
end
T1

T1 =

     0    40    40    40     0
    80    30    10    35   100
    80    20     0    25   100
    80    25     5    30   100
     0    20    20    20     0

% En funktion som beräknar nya temperaturer

% function Tn = pTemp(T)
%  [m,n] = size(T);
%  Tn = T;
%  for i = 2:m-1
%      for j = 2:n-1
%          Tn(i,j) = (T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1))/4;
%      end
%  end

T = T0;
T = pTemp(T)

T =

     0    40    40    40     0
    80    30    10    35   100
    80    20     0    25   100
    80    25     5    30   100
     0    20    20    20     0

% Vid villka temperaturer har systemet konvergerat?
T1 = T0; T = pTemp(T1);
while max(max(abs(T-T1)))>1e-3
    T1 = T;
    T = pTemp(T1);
end
T

T =

         0   40.0000   40.0000   40.0000         0
   80.0000   59.9986   56.7839   67.1415  100.0000
   80.0000   63.2125   59.9973   71.7839  100.0000
   80.0000   52.8558   48.2125   59.9986  100.0000
         0   20.0000   20.0000   20.0000         0

% jämför med svaret från lab2b
e = ones(9,1);
A = spdiags([-e -e 4*e -e -e],[-3 -1 0 1 3],9,9);
A(4,3) = 0; A(7,6) = 0; A(3,4) = 0; A(6,7) = 0;
b = [120; 40; 140; 80; 0; 100; 100; 20; 120];
x=A\b

x =

   60.0000
   56.7857
   67.1429
   63.2143
   60.0000
   71.7857
   52.8571
   48.2143
   60.0000

% imagesc
imagesc(T0)

imagesc(T0); colorbar

% Illustrera förloppet
T1 = T0; T = pTemp(T1);
while max(max(abs(T-T1)))>1e-3
    imagesc(T,[0 100]);
    T1 = T;
    T = pTemp(T1);
    pause(0.2);
end
disp('Klart');

Klart

% Använd finare grid
clear all; clc; clf;
n = 30;
T0 = zeros(n);
T0(2:n-1,1) = 80; T0(2:n-1,n) = 100;
T0(1,2:n-1) = 40; T0(n,2:n-1) = 20;

T1 = T0; T = pTemp(T1);
p = imagesc(T,[0,100]); colorbar;
while max(max(abs(T-T1)))>1e-1
    p.CData = T;  % uppdatera attribut istf att rita om hela figuren
    T1 = T;
    T = pTemp(T1);
    pause(0.1);
end
disp('Klart');

Klart

Ledning till Lab4

% Simulera en skogsbrand. Man tänker sig skogen indelad i små rutor:
% grön ruta -> kan börja brinna
% svart ruta -> kan ej brinna
% gul ruta -> brinner

% Modell för spridning:
% svart ruta -> förblir svart
% gul ruta -> blir svart i nästa tidssteg
% grön ruta -> blir gul med sannolikhet 0.9 om det finns en grannruta som �r
% gul

% Obs det finns ledning för hur man ska lösa lab4 (förberedande uppgifter).
% Så för att lösa labben:
% Läs igenom hela labpm, sätt dig in i de förberedande uppgifterna.

% Funktionen antalgrannar beräknar antal grannar som har värdet 2 i en
% matris

type antalgrannar

function antal = antalgrannar(A,i,j,n)
igrannar = max(1,i-1):min(n,i+1);
jgrannar = max(1,j-1):min(n,j+1);
if A(i,j) == 2
    antal = sum(sum(A(igrannar, jgrannar)==2))-1;
else
    antal = sum(sum(A(igrannar, jgrannar)==2));
end

clear all
A = [1 0 2;1 0 0;1 1 1]

A =

     1     0     2
     1     0     0
     1     1     1

antalgrannar(A,2,2,3)   % en brinnande granne

ans =

     1

antalgrannar(A,3,1,3)  % ingen brinnande granne

ans =

     0

Rita en startskog

figure(1); clf; hold on
n = 5;
for i = 1:n
    for j = 1:n
        F(i,j)=fill([i i+1 i+1 i],[j j j+1 j+1],'green');
    end
end

A = ones(n,n);  % Lagra alla gröna rutor

En svart och en gul ruta

A(1,1) = 0;
F(1,1).FaceColor='black';

A(3,4) = 2;
F(3,4).FaceColor='yellow';

% När vi är klara med hela programmet kommer vi lägga till
% axis off efter axis equal, så att skalorna på axlarna inte syns

Markera brinnande rutor

while 1
    [x,y,knapp] = ginput(1);
    if knapp ~= 1
        break;
    end
    x = floor(x); y = floor(y);
    F(x,y).FaceColor = 'yellow';
    A(x,y) = 2;
end

Simulera spridning av branden:

"loopa" igenom alla element i matrisen
Uppdatera till nästa tillstånd i en annan matris, B

% [m,n] = size(A);
% B = A;
% while finns2(A)
%     for i = 1:n
%         for j = 1:n
%             if A(i,j) == 2
%                 B(i,j) = 0;
%                 F(i,j).FaceColor='black';
%            ...
%         end
%     end
%     A = B;
%     pause(1);
% end

Published with MATLAB® R2016b