Bistable

Z Skrypty dla studentów Ekonofizyki UPGOW

(Różnice między wersjami)
(Utworzył nową stronę „===Obligacja ze stałym kuponem=== Mamy obligację, której emitent zobowiązuje się do płacenia odsetek regularnie raz do roku i zamierza zwrócić zaciągnięte …”)
m
 
(Nie pokazano 19 wersji pomiędzy niniejszymi.)
Linia 1: Linia 1:
-
===Obligacja ze stałym kuponem===
+
==Oscylator bistabliny==
-
Mamy obligację, której emitent zobowiązuje się do płacenia odsetek  regularnie raz do roku i zamierza zwrócić  zaciągnięte  zobowiązanie (wartość nominalną) w chwili wykupu, na koniec życia zobowiązania. Wartość takie obligacji dane jest [[IRF:Analiza_i_wycena_instrument%C3%B3w#Cena_godziwa_.28fair_price.29|wzorem]]
+
-
<math>\ P_o=\sum\limits_{i=1}^n\frac{C}{(1+r)^i} +\frac{P_N}{(1+r)^n},</math>
 
-
który możemy zaimplementować jako funkcję w matlabie:
 
-
<source lang="matlab">
+
===Układ ===
-
function P0=Bond_Fair_Price(PN,r,C,n)
+
-
  P0 = sum ( C./(1+r).^[1:n] ) + PN/(1+r)^n;
+
-
endfunction
+
-
</source>
+
-
[[Plik:fair_price.png|thumb|360px|Wykres zależności wartości obligacji od stopy procentowej]]
+
Rozważmy oscylator nieliniowy (usaa):
-
Dysponując taką funkcją możemy narysować wykres zależności ceny obligacji od stopy procentowej r:
+
-
<source lang="matlab">
+
-
r=linspace(0.01,0.5,25)
+
-
for i=1:length(r)
+
-
  p(i)=Bond_Fair_Price(100,r(i),7,25);
+
-
endfor
+
-
plot(r,p)
+
-
</source>
+
-
''Uwaga: Proszę zwrócić uwagę na fragment:''
+
<math>m \ddot x = -\frac{dU(x)}{dx}-\gamma x  ,</math>  
-
<source lang="matlab">
+
-
  C./(1+r).^[1:(n-1)]
+
-
</source>
+
-
''który tworzy wektor o elementach będących funkcją wskaźnika''
+
-
<math>\frac{C}{(1+r)^i} </math> dla <math>i=1..(n-1)</math>.
+
równoważnie:
-
''Działanie można sobie przećwiczyć na prostym przykładzie:''
+
<math>\dot x = v</math>  
-
<source lang="matlab">
+
-
octave:271> idx=[1:4]
+
-
idx =
+
-
  1  2  3  4
+
<math>m \dot v = -\frac{dU(x)}{dx}-\gamma x  ,</math>
-
octave:272> idx.^2
+
jako potencjal weżmy funkcję z dwoma minimami np.:
-
ans =
+
-
    1    4   9  16
+
<math>U(x) = \displaystyle x^4-4 x^2,</math>
 +
 
 +
Wykres U(x) można otrzymać poleceniem:
 +
<source lang="matlab">
 +
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>
 +
<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);
-
Przeliczmy teraz [[IRF:Analiza_i_wycena_instrument%C3%B3w#Cena_godziwa_.28fair_price.29|przykład]] ze skryptu Instrumenty Rynku :
+
  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===
<source lang="matlab">
<source lang="matlab">
-
octave:157>P0=1
+
lsode_options("absolute tolerance",1e-6)
-
octave:157>PN=100
+
lsode_options("relative tolerance",1e-6)
-
octave:157>r=0.07
+
n=500;
-
octave:157>C=6
+
tmax=10;
-
octave:157>Bond_Fair_Price(PN,r,C,2)
+
t = linspace(0,tmax,n);
-
ans =  98.192
+
global gama;
 +
gama=0.9
 +
ny=10;
 +
y=linspace(0,5,ny);
 +
x0=-sqrt(2);
 +
trajsxP=zeros(ny,n);
 +
trajsyP=zeros(ny,n);
 +
trajsxN=zeros(ny,n);
 +
trajsyN=zeros(ny,n);
 +
for iy=1:ny
 +
  X0 = [x0,y(iy)];
 +
  [X,T,MSG]=lsode(@ODEbistable,X0,t);
 +
  if(X(length(X),1)>0)
 +
    trajsxP(iy,:)=X(:,1);
 +
    trajsyP(iy,:)=X(:,2);
 +
  else
 +
    trajsxN(iy,:)=X(:,1);
 +
    trajsyN(iy,:)=X(:,2);
 +
  endif
 +
endfor
 +
disp("Integration done");
 +
hold on
 +
every=20
 +
for i=(every+1):every:size(trajsxP)(2)
 +
  plot(trajsxP(:,(i-every):i)',trajsyP(:,(i-every):i)',"r-",\
 +
    trajsxN(:,(i-every):i)',trajsyN(:,(i-every):i)',"b-")
 +
  drawnow()
 +
  disp(t(i));
 +
  fflush(stdout);
 +
  print ( sprintf("Traj_anim_%05d.png",i), "-S300,220");
 +
endfor
 +
hold off
</source>
</source>
-
Analogicznie możemy zaimplementować naszą funkcje dla przypadków: [[IRF:Analiza_i_wycena_instrument%C3%B3w#Obligacja__zerokuponowa|m wypłat kuponu]] w jednym roku:
+
który możemy zaimplementować jako funkcję w matlabie:
-
<source lang="matlab">
+
-
function P0=Bond_Fair_Price_multi(PN,r,C,n,m)
+
-
  P0 = sum ( (C/m)./(1+r/m).^[1:n] ) + PN/(1+r/m)^n;
+
-
endfunction
+
-
</source>
+
-
oraz dla  [[IRF:Analiza_i_wycena_instrument%C3%B3w#Wycena_przy_kapitalizacji_ci.C4.85g.C5.82ej|kapitalizacji ciągłej]] gdzie mamy:
 
<source lang="matlab">
<source lang="matlab">
-
function P0=Bond_Fair_Price_cont(PN,r,C,t)
+
function dx = ODEbistable(X,T)
-
ti=1:floor(t);
+
    global gama;
-
P0 =  sum (C*exp(-r*ti)) + (t-floor(t))*C*exp(-r*t) + PN*exp(-r*t);
+
    dx = zeros(2,1);
-
endfunction
+
    dx(1) = X(2);
 +
    dx(2) = -gama*X(2)-4*X(1).^3+8*X(1);
 +
    return
 +
end
</source>
</source>
-
[[Plik:fair_price_cont_discr.png|thumb|360px|Wykres zależności wartości obligacji od stopy procentowej: porównanie modeli z kapitalizacjami: ciągłą i dyskretną.]]
+
[[Plik:bistable_anim.gif|thumb|216px|Trajektorie w nieliniowym oscylatorze.]]
 +
[[Plik:bistable_color.png|thumb|360px|Baseny przyciągania w nieliniowym oscylatorze.]]

Aktualna wersja na dzień 09:00, 8 lis 2010

Oscylator bistabliny

Układ

Rozważmy oscylator nieliniowy (usaa):

\(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

lsode_options("absolute tolerance",1e-6)
lsode_options("relative tolerance",1e-6)
n=500;
tmax=10;
t = linspace(0,tmax,n); 
global gama;
gama=0.9
ny=10;
y=linspace(0,5,ny); 
x0=-sqrt(2);
trajsxP=zeros(ny,n);
trajsyP=zeros(ny,n);
trajsxN=zeros(ny,n);
trajsyN=zeros(ny,n);
for iy=1:ny
  X0 = [x0,y(iy)];
  [X,T,MSG]=lsode(@ODEbistable,X0,t);
  if(X(length(X),1)>0)
    trajsxP(iy,:)=X(:,1);
    trajsyP(iy,:)=X(:,2);
  else
    trajsxN(iy,:)=X(:,1);
    trajsyN(iy,:)=X(:,2);
  endif
endfor
disp("Integration done");
hold on
every=20
for i=(every+1):every:size(trajsxP)(2)
  plot(trajsxP(:,(i-every):i)',trajsyP(:,(i-every):i)',"r-",\
     trajsxN(:,(i-every):i)',trajsyN(:,(i-every):i)',"b-")
  drawnow()
  disp(t(i));
  fflush(stdout);
  print ( sprintf("Traj_anim_%05d.png",i), "-S300,220");
endfor
hold off

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
Trajektorie w nieliniowym oscylatorze.
Baseny przyciągania w nieliniowym oscylatorze.
Źródło „http://172.16.3.1/ekonofizyka/index.php/Bistable”