Guest User

Untitled

a guest
Sep 15th, 2019
96
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