Bistable
Z Skrypty dla studentów Ekonofizyki UPGOW
(Różnice między wersjami)
(→Układ) |
(→Układ) |
||
Linia 20: | Linia 20: | ||
fplot(@(x) x.^4-4*x.^2,[-2.1,2.1],200) | fplot(@(x) x.^4-4*x.^2,[-2.1,2.1],200) | ||
</source> | </source> | ||
- | |||
<math>f(x) =-\frac{dU(x)}{dx} = -4 x^3+8 x^2,</math> | <math>f(x) =-\frac{dU(x)}{dx} = -4 x^3+8 x^2,</math> | ||
+ | <source lang="matlab"> | ||
+ | close all | ||
+ | clear all | ||
+ | f =@(x,E) sqrt(8*x.^2-2*x.^4-E); | ||
+ | xa=linspace(-2,2,200); | ||
+ | xb=linspace(2,-2,200); | ||
+ | hold on | ||
+ | for E=-8:2:8 | ||
+ | Rts=roots([-2,0,8,0,-E]); | ||
+ | AllRealRoots=sort(Rts(find (imag(Rts)==0))); | ||
+ | if (length(AllRealRoots)==2) | ||
+ | xminmax=sort([AllRealRoots',-0.00001,0.00001]); | ||
+ | else | ||
+ | xminmax=AllRealRoots'; | ||
+ | endif | ||
+ | xa=linspace(xminmax(1),xminmax(2),200); | ||
+ | xb=linspace(xminmax(3),xminmax(4),200); | ||
+ | |||
+ | if (E==0) | ||
+ | plot([rotdim(xa,2),xa],[-f(rotdim(xa,2),E),f(xa,E)],"r-;E=0;"); | ||
+ | plot([xb,rotdim(xb,2)],[f(xb,E),-f(rotdim(xb,2),E)],"r-"); | ||
+ | else | ||
+ | plot([rotdim(xa,2),xa],[-f(rotdim(xa,2),E),f(xa,E)],"g-"); | ||
+ | plot([xb,rotdim(xb,2)],[f(xb,E),-f(rotdim(xb,2),E)],"g-"); | ||
+ | endif | ||
+ | endfor | ||
+ | hold off | ||
+ | </source> | ||
===Analiza=== | ===Analiza=== |
Wersja z 10:16, 27 paź 2010
Oscylator bistabliny
Układ
Rozważmy oscylator nieliniowy:
\(m \ddot x = -\frac{dU(x)}{dx}-\gamma x ,\)
równoważnie:
\(\dot x = v\)
\(m \dot v = -\frac{dU(x)}{dx}-\gamma x ,\)
jako potencjal weżmy funkcję z dwoma minimami np.:
\(U(x) = \displaystyle x^4-4 x^2,\)
Wykres U(x) można otrzymać poleceniem:
fplot(@(x) x.^4-4*x.^2,[-2.1,2.1],200)
\(f(x) =-\frac{dU(x)}{dx} = -4 x^3+8 x^2,\)
close all clear all f =@(x,E) sqrt(8*x.^2-2*x.^4-E); xa=linspace(-2,2,200); xb=linspace(2,-2,200); hold on for E=-8:2:8 Rts=roots([-2,0,8,0,-E]); AllRealRoots=sort(Rts(find (imag(Rts)==0))); if (length(AllRealRoots)==2) xminmax=sort([AllRealRoots',-0.00001,0.00001]); else xminmax=AllRealRoots'; endif xa=linspace(xminmax(1),xminmax(2),200); xb=linspace(xminmax(3),xminmax(4),200); if (E==0) plot([rotdim(xa,2),xa],[-f(rotdim(xa,2),E),f(xa,E)],"r-;E=0;"); plot([xb,rotdim(xb,2)],[f(xb,E),-f(rotdim(xb,2),E)],"r-"); else plot([rotdim(xa,2),xa],[-f(rotdim(xa,2),E),f(xa,E)],"g-"); plot([xb,rotdim(xb,2)],[f(xb,E),-f(rotdim(xb,2),E)],"g-"); endif endfor hold off
Analiza
który możemy zaimplementować jako funkcję w matlabie:
function dx = ODEbistable(X,T) global gama; dx = zeros(2,1); dx(1) = X(2); dx(2) = -gama*X(2)-4*X(1).^3+8*X(1); return end