Równania różniczkowe z laboratorium - wybrane kody przykładowe

Poniżej znajdą Państwo kody źródłowe podsumowujące tematy omawiane przez nas na ćwiczeniach w laboratorium. Nie jest to kopia tego, co robiliśmy na ćwiczeniach. Nie jest to też stricte materiał dydaktyczny - raczej coś w rodzaju "bloga z kodami", ot dla zgrubnej informacji, co było robione.

Poniższe pliki źródłowe można pobierać, klikając w nazwę programu.

Warto zmienic nazwe pliku tak, by pozbyc sie przedrostka, czyli np. zamiast c5-harmosc.m zapisac pod nazwa harmosc.m

Zestaw 6

c6-impliciteuler.m

  1. span style="color: #228B22;">% Dla równania $x' = f(t,x), \quad t_0 < t \leq T$ z warunkiem początkowym $x(t_0) = x_0$  
  2. % wyznacza x(T)
  3. % metodą niejawną Eulera z krokiem $h = (T-t_0)/N$
  4. % nowy x to x(t+h)

c6-stiffness.m

  1. span style="color: #228B22;">% jesteśmy więc b. blisko punktu stacjonarnego, 
  2. % gdzie wartości własne to 0 (skąd degeneracja Jakobianu i problem w fsolve) 
  3. % oraz dwie ujemne: -0.405 oraz -2189.6
  4. % ... dla Eulera
  5. % ... dla iEulera
  6. % wykonujemy JEDEN krok schematu...

Zestaw 5

c5-harmosc.m

  1. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  2. %
  3. % oscylator harmoniczny x'' + 2*p*x' + w0^2*x = F(t)
  4. %
  5. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  6.  
  7.  
  8. p  = 0;		% tłumienie
  9. w0 = 6;		% częstość własna
  10. A  = 3;		% amplituda wymuszenia
  11. w  = 4;		% częstość wymuszenia
  12. T  = 100;	% czas obserwacji
  13. % zewnętrzne wymuszenie
  14.  
  15. % zamieniamy na układ równań:
  16. %
  17. % x' = v
  18. % v' = F(t) - 2*p*v - w0^2*x
  19. %
  20. % X = [X(1), X(2)] = [x,v]
  21. % wykres rozwiązania
  22. % wykres fazowy

Zestaw 4

c4-euler.m

  1. % przed użyciem funkcji należy ją zapisać do pliku o nazwie `euler.m`
  2.  
  3. %%%%%%%%%%%%%%%%
  4. %
  5. % Metoda Eulera 
  6. %
  7. %%%%%%%%%%%%%%%%
  8. % Dla równania $x' = f(t,x), \quad t_0 < t \leq T$ z warunkiem początkowym $x(t_0) = x_0$  
  9. % wyznacza x(T)
  10. % metodą Eulera z krokiem $h = (T-t_0)/N$
  11. % nowy x to x(t+h)

c4-plotsolution.m

  1. %%%%%%%%%%%%%%%%%%%%%%
  2. %
  3. % Testowanie schematów 
  4. %
  5. %%%%%%%%%%%%%%%%%%%%%%
  6.  
  7.  
  8. % definicja problemu 
  9. %
  10. % prawa strona równania, $f(t,x)$
  11. t0 = 0; x0 = 1; T = 4; % dane zadania
  12.  
  13. % definicja rozwiązania dokładnego, $x(T)$ - dla sprawdzenia błędu
  14. %
  15. % definicja parametrów schematu
  16. %
  17. N = 20; % liczba kroków 
  18.  
  19. %
  20. % cz. I: test schematu
  21. %
  22. %g\n\t Euler: %f | RK4: %f | dokł: %f\n', T, xe, xrk4, X(T)); % wartości w chwili T
  23. %g\n\t Euler: %e | RK4: %e | dokł: %e\n', T, [xe, xrk4,X(T)]-X(T)); % błędy w chwili T
  24.  
  25. %
  26. % cz. II: wykres rozwiązania
  27. %
  28. % punkty do wykresu
  29. % wartości rozwiązania w chwilach $t_0,\ldots, T$ ...
  30. % ... dla Eulera
  31. xrk4 = xe; 							% ... dla RK4
  32. % wykonujemy JEDEN krok schematu...
  33. %g\n\t Euler: %f | RK4: %f | dokł: %f\n', T, xe(end), xrk4(end), X(T)); % wartości w chwili T
  34. %g\n\t Euler: %e | RK4: %e | dokł: %e\n', T, [xe(end), xrk4(end),X(T)]-X(T)); % błędy w chwili T
  35.  
  36. % rysujemy wykresy
  37. %
  38. % cz. III: sprawdzenie rzędu
  39. %
  40. % uwaga: jeśli weźmiemy za dużo, może okazać się, że rząd 'nie wychodzi'
  41. %f', log2(err_e(1:end-1)./err_e(2:end)));
  42. %f', log2(err_rk4(1:end-1)./err_rk4(2:end)));

c4-rk4.m

  1. % przed użyciem funkcji należy ją zapisać do pliku o nazwie `rk4.m`
  2.  
  3. %%%%%%%%%%%%%%%%
  4. %
  5. % Metoda RK4 
  6. %
  7. %%%%%%%%%%%%%%%%
  8. % Dla równania $x' = f(t,x), \quad t_0 < t \leq T$ z warunkiem początkowym $x(t_0) = x_0$  
  9. % wyznacza x(T)
  10. % metodą RK4 z krokiem $h = (T-t_0)/N$
  11. % nowy x to x(t+h)

Zestaw 3

c3-maxima-ode-numer.mac

  1. span style="color: #808080; font-style: italic;">;
  2. ; 		/* definiujemy równanie */
  3. ;			/* wersja minimalistyczna */
  4. ; 	/* równanie z dwoma parametrami - czy można je rozwiązać symbolicznie? */
  5. "k=1,a=1""k=-1:1,a=-3:3"]);
  6.  

Zestaw 2

c2-maxima-ode-symbol.mac

  1. span style="color: #808080; font-style: italic;">; 	/* definiujemy równanie */
  2. ; 	/* rozwiązujemy */
  3. method; 				/* sprawdzamy typ-metodę */
  4. ;		/* wyznaczamy stałe dla zadania Cauchy'ego */ 
  5. ;	/* sprawdzenie: podstawiamy do równania rozwiązanie... */
  6. ;		/*    ...przykładamy pochodną... */
  7. ;
  8. ; /* seria rozwiązań ze zmienną stałą %c */
  9. ; 
  10. ;	/* rozwiązanie, wynik w  postaci uwikłanej */
  11. depends(x,t); 				/* koniecznie informujemy, że x jest funkcją t! */
  12. ; 			/* różniczkujemy uwikłane wyrażenie */
  13. ;		/* rozwikłujemy względem pochodnej i ewentualnie upraszczamy */
  14. ;					/* zapominamy, czym jest x */
  15. ;
  16. ; /* cała seria wartości stałych */
  17. ;			/* pakiet zawierający więcej metod rozwiązywania */
  18. put('contrib_ode,true,'verbose); /* więcej diagnostyki */
  19. f(dx) := dx^2;
  20. ; /* definiujemy równanie */
  21. ; /*...i rozwiązanie, które tym razem jest listą */
  22. method;
  23. ;		/* ...i które łatwiej sprawdzić */ 
  24. ;	/* uwzględniamy warunek początkowy */
  25. ;		/* czyści pamięć */
  26.  

Zestaw 1

c1-maxima-intro.mac

  1. span style="color: #808080; font-style: italic;">;
  2. ; /* arccos(1/2) nie zadziała! */
  3. ;
  4. ;
  5. ;
  6. ;
  7. ;
  8. ;			/* wartość numeryczna */
  9. ;	/* sprawdzamy rozwiązanie: podstawiamy do równania... */
  10. ;					/* ...rozwijamy */
  11. ;
  12. ;
  13. da: ratsimp(da);
  14.  
  15. /* całkowanie */
  16. c: integrate(x^4/(x+n),x); 		/* nieoznaczona! */
  17. c: integrate(x^4/(x+n),x,0,%pi); /* oznaczona! */
  18. ;
  19. c: c, n=4;  /* to samo inaczej i na skróty */
  20. ;
  21. ; 			/* obcinamy za duże wartości */
  22. ;	/* dwa wykresy na jednym */
  23. ;		/* czyści pamięć */
  24.  
Piotr Krzyżanowski | Aktualizacja: 27.02.2019, 10:39:56.