SHARE
TWEET

Untitled

a guest Sep 15th, 2019 83 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. % Laboration i FMM
  2. %% Uppgift 1 & 2
  3. clear all
  4. % close all
  5. load randwebs
  6. load randinits
  7.  
  8. IVP = 225;        %vilka randweb & randinits som ska köras
  9. A = randwebscell{IVP};
  10. x0 = randinitscell{IVP};
  11. NoNeg = odeset('NonNegative',1);
  12.  
  13. tmax = 1500;
  14. tspan = [1 tmax];
  15.  
  16. [t, x] = (ode45(@(t,x) x.*A*x,tspan,x0,NoNeg)) ;
  17.  
  18. plot(t, x)
  19.  
  20.  
  21. n = 0;
  22. for i = 1:10
  23.     if x(end, i) > 10^(-4)
  24.         n= n+1;
  25.     end
  26. end
  27. disp('Number of species alive')
  28. n
  29. %% Uppgift 3
  30. alive = zeros(1,1000);
  31. for i = 1:1000
  32.     i
  33.     A = randwebscell{i};
  34.     X0 = randinitscell{i};
  35.     [t,x] = ode45(@(t,x)x.*A*x,tspan,X0,odeset('NonNegative',1));
  36.     for j = 1:10
  37.         if x(end,j) > 10^-3
  38.             alive(i) = alive(i) + 1;
  39.         end
  40.     end
  41. end
  42. %% Uppgift 4
  43. histogram(alive)
  44.  
  45.  
  46. %% uppgift 9 (att oscillera eller inte oscillera, det är frågan)
  47.     % Den här testar bara en i taget, det var så vi jobbade när du var här
  48.     % Se nästa för färdig version
  49.  
  50.  
  51. grafNr = 10;
  52. annat = A(grafNr, :)*x';
  53. dx = (x(:,grafNr)).*annat';
  54. figure
  55. plot(t, dx)
  56.  
  57. if isempty(find(abs(dx(round((0.75*length(dx))):length(dx))) > 10^(-4))) == 1
  58.     disp('Skiten oscillerar inte')
  59. else
  60.     disp('Oscillerar')
  61. end
  62.  
  63.  
  64. %% Riktiga uppgift 9, räkna hur många som oscillerar
  65.     %verkar funka, men dubbelkolla gärna
  66.  
  67. n = 0;
  68. for i = 1:10
  69.     grafNr = i;
  70.     annat = A(grafNr, :)*x';
  71.     dx = (x(:,grafNr)).*annat';
  72.    
  73.    if isempty(find(abs(dx(round(0.75*(length(dx)))):length(dx))) > 10^(-4)) == 0
  74.        n=n+1;
  75.    end
  76. end
  77.  
  78. disp('Antal oscillerande: ')
  79. n
  80.  
  81.  
  82. %% Uppgift 3 & 4
  83. % Max tidsintervall?
  84. tspan = [1 2500];
  85.  
  86.  
  87. alive = zeros(1,1000);
  88. sensitivity = 10^(-3)
  89.  
  90. for j= 1:1000
  91.     A = randwebscell{j};
  92.     x0 = randinitscell{j};
  93.    
  94.     [t, x] = (ode45(@(t,x) x.*A*x,tspan,x0,NoNeg)) ;
  95.    
  96.     for i = 1:10
  97.         if x(end, i) > sensitivity
  98.             alive(j) = alive(j) + 1;
  99.         end
  100.     end
  101. end
  102.  
  103. figure
  104. histogram(alive)
  105.  
  106. %Hur många har 10 levande?
  107. tens = find(alive == 10)
  108. nrOf10 = length(tens)
  109.  
  110. % Uppgift 5, plotta 10 överlevare
  111. A = randwebscell{225};
  112. x0 = randinitscell{225};
  113. [t, x] = (ode45(@(t,x) x.*A*x,tspan,x0,NoNeg)) ;
  114. figure
  115. plot(t, x)
  116. hold on
  117. title('Nr 225, två som återuppstår vid t = 1500')
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top