ریاضیات برای اقتصاد/بهینه‌سازی پویا

از testwiki
پرش به ناوبری پرش به جستجو

الگو:سرصفحه تا به اینجای کتاب صرفاً به بررسی مباحث در یک مقطع از زمان پرداختیم. سوالاتی که به آن‌ها پاسخ دادیم حداکثرسازی مطلوبیت فرد، یا بنگاه تنها در یک مقطع از زمان بود. در طول زمان ممکن است قیمت‌ها، کالاهای مصرفی افراد، نهاده‌های مورد مصرف بنگاه‌ها، تابع مطلوبیت افراد یا تابع تولید بنگاه‌ها تغییر کنند. در این بخش بیشتر به سوالات این‌چنینی می‌پردازیم و به‌جای یک نقطه در زمان، در یک بازه از زمان و یک افق زمانی اقدام به حداکثر یا حداقل‌سازی می‌کنیم.

هرگاه بحث افق‌های زمانی شود، از معادلات دیفرانسیل استفاده می‌کنیم. تغییر در زمان، علامت تغییر، و خود استخراج خود متغیر با استفاده از معادلات دیفرانسیل شدنی است. در حال حاضر ما صرفاً با معادلات دیفرانسیل معمولی الگو:انگلیسی (فقط با یک متغیر و آن هم‌زمان) سر و کار داریم و نه با الگو:انگلیسی که توابع چند متغیره را بررسی می‌کند.

هدف ما از پیگیری این موضوع بیشتر رسم نگاره‌ٔ فاز، یا صفحهٔ فاز الگو:انگلیسی است. نگارهٔ فاز، الگو:انگلیسی کل میدان برداری است و مسیر فاز الگو:انگلیسی، خطوطی است که مسیر سیستم با گذر زمان را معین می‌کنند.[۱]

پیش مطالعه

یادآوری

توابع نمایی و مشتق آن‌ها

الگو:تمرین هر تابع نمایی را می‌توان با پایهٔ عدد اویلر (یا هر عدد دیگری) بیان کرد. الگو:چپ‌چین 2t=(eln(2))t=eln(2)t الگو:پایان از ویدئوی یادآوری بالا به خاطر داریم که مشتق هر تابع نمایی خود آن تابع ضرب در یک عدد ثابت است، که آن عدد ثابت با لگاریتم طبیعی تعیین می‌شود؛ داریم:

ddt2t=ddteln(2)t=ln(2)eln(2)t=ln(2)2t الگو:پایان تمرین

تعاریف

یک سیستم از معادلات دیفرانسیل معمولی و خطی به این صورت نشان داده می‌شود: x˙=Ax که در آن 𝐱˙=d𝐱dt=[x˙1x˙n]=[dx1dtdxndt]

جواب چنین دستگاهی برابر با 𝐱(t)=eAt𝐱0 است که در آن 𝐱0 نقطهٔ آغاز دستگاه است و خود متغیر بر حسب زمان نوشته شده است. چنین جوابی از آنجا به دست آمد که تنها تابعی که مشتقش برابر با خودش در یک ضریب است، همان تابع نمایی e است. الگو:تمرین n=1xx˙=axx(t)=eatx0

که در آن x0 یک ثابت است و نقطهٔ شروع نیز هست زیرا به ازای t=0 عبارت x(t)=eatx0 برابر با x0 می‌شود. الگو:پایان تمرین الگو:تمرین نکته: ex=1+x+x22!+x33!++xnn!+=i=01n!xn[۲] دیدن این فیلم برای درک بهتر این بسط توصیه می‌شود. الگو:پایان تمرین

  • دستگاه غیرمزدوج الگو:انگلیسی: مشتق هر متغیر تابعی بر حسب خود آن متغیر است و نه سایر متغیرها. مثلاً x˙1=1x1 و x˙2=2x2 معادل ماتریسی این گونه دستگاه‌ها قطری می‌شود.

مثال: [x˙1x˙2]=[1002][x1x2]

  • سامانه پویا الگو:انگلیسی یک مدل‌سازی ریاضی برای هرگونه قانون ثابت است که وابستگی موقعیت یک نقطه در یک فضای پیراگیر نسبت به یک پارامتر (که معمولا زمان است) را توضیح می‌دهد.[۳] سامانه پویا ابزاری است که به ما کمک می‌کند تغییر وضعیت از یک حالت به حالت دیگر را بررسی کنیم.[۴] سامانه پویا می‌تواند به صورت یک تابع از ×2 (یک متغیرt و دو متغیر نمایانگر نقطه‌ی شروع همچون [c1c2] به عنوان دامنه) به 2 (مکان نهایی نقطه‌ی اولیه‌ی c در زمان t) تعریف شود.

سامانه‌ی پویای مثال بالا چنین است: ϕ(t,c)=[et00e2t]𝐜 که همانطور که در بالا گفته شد، در آن c یک ماتریس و نمایانگر نقطه‌ی آغاز است، حاصل ضرب نیز (مقدار تابع) مکان نهایی نقطه‌ی آغاز خواهد بود. ضمنا همانطور که مشاهده می‌شود برای دستیابی به سامانه‌ی پویا نیاز به حل دستگاه داریم، اما حل دستگاه معادلات دیفرانسیل در اغلب اوقات دشوار است. از این رو می‌توان به‌جای سامانه‌ی پویا و بدون حل دستگاه، با خود دستگاه اولیه‌ای که داشتیم اقدام به رسم f(x)=Ax کنیم که نگاشتی با تعریف f:2×2 است. دامنه‌ی این نگاشت هر x2 و برد آن بردارهایی همچون f(x) است که به ما جهت و اندازه‌ی حرکت در نقطه‌ی x به ازای تغییرات t را می‌دهد.

رسم نگارهٔ فاز: مباحث اولیه

رسم نگارهٔ فاز با گنو اکتاو

نگارهٔ فاز ساده‌ی رسم شده برای مثال معرفی‌شده، با استفاده از گنو اکتاو

برای رسم نمودار فاز دستگاه غیر مزدوج

[x˙1x˙2]=[1002][x1x2]

، علاوه بر راهکارهای ریاضیاتی که در پیش رو معرفی می‌شود می‌توان از نرم‌افزار GNU Octave[۵] استفاده کرد:

figure;

[x1,x2]=meshgrid(-2:0.2:2,-2:0.2:2);

dx1 = -x1;		%x1'= -x
dx2 = x2.*2;		%x2'= 2x2

quiver(x1,x2,dx1,dx2) 	%plot the vector field with arrows
xlabel('x1');
ylabel('x2');
title("Phase Plot for The Wikibooks Example");
xlim([-2 2]);
ylim([-2 2]);

grid on;
hold on;

رسم نگارهٔ فاز با استفاده از ریاضیات صرف

نمودار فوق چگونه رسم شده است؟ دو نقطهٔ (1,1) و (0,1) را در نظر بگیرید. به ازای این دو نقطه داریم: [1002][01]=[02][1002][11]=[12] همانطور که دیده می‌شود، صرفا با داشتن دستگاه و بدون حل آن توانستیم بردارهای مماس بر مسیر فاز را رسم کنیم. البته اندازهٔ این بردارها همگی با عددی کوچک‌تر مقیاس می‌شود تا جا برای رسم سایر بردارها نیز باقی بماند.

رسم مسیر فاز در کنار نگارهٔ فاز

برای اینکه مسیر فاز را نیز داشته باشیم داریم:[۶]

f=@(t,X) [-1*X(1);2*X(2)];
a=[1,-0.1,2,-2];
b=[-0.2,0.1,0.1,-0.1];
for k=1:4
  [t,Sap]=ode45(f,[0 2],[a(k),b(k)]);
  plot(Sap(:,1),Sap(:,2),'LineWidth',2,'color','k')
end

که افزودن این قطعه کد به کد پیشین، به ما چنین نموداری خواهد داد:

مسیر فاز به‌ازای خوراندن نقاط اولیهٔ (1,0.2)،(0.1,0.1)،(2,0.1)،(2,0.1) به دستگاه، در کنار فضای برداری (نگارهٔ فاز)

در این قسمت از درس می‌توان بیش از پیش به اهمیت بردارهای ویژه پی برد. از آنجایی که ماتریس اولیهٔ ما قطری بود، بردارهای ویژه‌اش همان پایه‌های استاندارد هستند (یعنی بردارهای ویژه هم‌جهت با محورهای

x1

و

x2

هستند) همانطور که می‌بینید هر نقطه‌ای روی این دو محور صرفا مقیاس الگو:انگلیسی شده و تغییر جهت نمی‌دهد.

بررسی بهتر اهمیت بردارهای ویژه

بیاید با هم بردارهای سازنده‌ی نگارهٔ فاز بالا را با ماتریس دوران، ۴۵ درجه چرخش دهیم، به‌گونه‌ای که ماتریس، مقادیر ویژه‌اش تغییر نکند ولی بردارهای ویژه‌اش ۴۵ درجه تغییر زاویه دهند (یک ماتریس مشابه ماتریس [1002] را با استفاده از ماتریس دوران ۴۵ درجه‌ای و معکوسش می‌سازیم.)

B=P1AP=[12121212][1002][1111]=[121121][1111]=[12323212]

دستگاه معادلات دیفرانسیل تبدیل یافته چنین خواهد بود: x˙=Bx یا به عبارتی: [x˙1x˙2]=[12323212][x1x2]

محاسبهٔ مقادیر و بردارهای ویژه با گنو اکتاو

چه از طریق حل با قلم و کاغذ و چه با اکتاو (پس از نصب پکیج Symbolic)، می‌توان برابری مقادیر ویژه با ماتریس قطری اولیه را آزمود:

pkg load symb;
A=sym([0.5,-1.5;-1.5,0.5]);
[V,D]=eig(A)

جواب این کد چنین است:

V = (sym 2×2 matrix)

  ⎡1  -1⎤
  ⎣1  1 ⎦

D = (sym 2×2 matrix)

  ⎡-1  0⎤
  ⎣0   2⎦

مشاهده می‌شود که مقادیر ویژه تغییر نیافته‌اند، ولی بردار ویژه‌ی متناظر با مقدار ویژه‌ی

1

از

e1=[10]

به

Pe2=[1111][10]=[11]

تبدیل یافته است. بردار ویژه‌ی متناظر با

2

نیز، از

e2=[01]

به

Pe2=[1111][01]=[11]

تبدیل یافته است.

درک اهمیت بردارهای ویژه: رسم مسیر بردارهای ویژه در کنار نگارهٔ فاز با گنو اکتاو

مشابه قبل، نگارهٔ فاز دستگاه جدید که ۴۵ درجه چرخیده است، یعنی [x˙1x˙2]=[12323212][x1x2] را رسم می‌کنیم:

برای دیدن کد گسترش را بزنید.

[x1,x2]=meshgrid(-2:0.2:2,-2:0.2:2);
figure;
dx1= 0.5.*x1-1.5.*x2;
dx2= -1.5.*x1+0.5.*x2;
quiver(x1,x2,dx1,dx2);
hold on;
f=@(t,X) [1/2*X(1)+-3/2*X(2);-3/2*X(1)+1/2*X(2)];
a=[1.8,-2,0,-0.3];
b=[1.5,-1.8,0.2,-0.5];
for k=1:4
  [t,Sap]=ode45(f,[0 2],[a(k),b(k)]);
  plot(Sap(:,1),Sap(:,2),'LineWidth',2,'color','k')
end

g=-2:0.2:2;
plot(g,g,'r');
plot(g,-g,'r');
axis equal;
xlabel('x1');
ylabel('x2');
title("Phase portrait of rotated system of equation of the Wikibooks example");
xlim([-2 2]);
ylim([-2 2]);
نگارهٔ فاز دستگاه معادلات دیفرانسیل جدید که ۴۵ درجه چرخیده است. دو خط قرمز هم‌جهت با دو بردار ویژه هستند.

همینطور برای مسیر فاز این بردار رسم زیر را داریم:

clear
[x1,x2]=meshgrid(-2:0.2:2,-2:0.2:2);
figure;
dx1= 0.5.*x1-1.5.*x2;
dx2= -1.5.*x1+0.5.*x2;
quiver(x1,x2,dx1,dx2);
hold on;
f=@(t,X) [1/2*X(1)+-3/2*X(2);-3/2*X(1)+1/2*X(2)];
a=[1.8,-2,0,-0.3];
b=[1.5,-1.8,0.2,-0.5];
for k=1:4
  [t,Sap]=ode45(f,[0 2],[a(k),b(k)]);
  plot(Sap(:,1),Sap(:,2),'LineWidth',2,'color','k')
end

g=-2:0.2:2;
plot(g,g,'r');
plot(g,-g,'r');
axis equal;
xlabel('x1');
ylabel('x2');
title("Phase portrait of rotated system of equation of the Wikibooks example");
xlim([-2 2]);
ylim([-2 2]);

grid on;
نمودار فاز چرخیده به همراه چهار مسیر فاز در نقاط (1.8,1.5) و (2,1.8) و (0,0.2) و (0.3,0.5). اهمیت بردارهای ویژه (که بر روی مسیر خطوط قرمز رنگ قرار می‌گیرند) در اینجا بهتر از پیش درک می‌شود.

حل دستگاه‌های مزدوج

فرض کنید که x˙=Ax دستگاهی مزدوج باشد. در این صورت حل آن سخت است. با استفاده از قطری‌سازی می‌توان دستگاه مزدوج را غیر مزدوج ساخته و راحت‌تر به جواب آن دست یافت.

اثبات راهکار

ابتدا ماتریس P=[𝐯1𝐯2𝐯n] را با استفاده از بردارهای ویژه‌ی A تشکیل می‌دهیم. سپس وارون این ماتریس یعنی P1 را به دست می‌آوریم. متغیری همچون y=P1x را تعریف می‌کنیم. با گرفتن مشتق هر دو طرف داریم: y˙=P1x˙ از طرفی خود دستگاه اولیه هست: x˙=Ax پس با جایگذاری x˙=Ax داریم: y˙=P1Ax از نتیجه‌ی تعریف متغیر y=P1x می‌توان Py=x را نتیجه گرفت، با جایگذاری این نتیجه در معادلهٔ قبلی داریم: y˙=P1APy. دستگاه به دست آمده قطری است و مقادیر روی قطر آن مقادیر ویژه‌ی ماتریس A هستند یعنی: y˙=diag[λ1,,λn]y. چنین دستگاهی دیگر مزدوج نیست و به عنوان یک دستگاه غیرمزدوج راه حل کلی y(t)=diag[eλ1t,,eλnt]y(0) را دارد.

داشتیم: Py=x، از این معادله ما y را به دست آورده‌ایم. می‌توان دو طرف تساوی y(t)=diag[eλ1t,,eλnt]y(0) را از چپ در P ضرب کرد در این صورت داریم: Py(t)=x(t)=Pdiag[eλ1t,,eλnt]y(0)، از طرفی طبق تعریف خود y=P1x می‌توان y(0)=P1x(0) را نیز داشت اگر کل diag[eλ1t,,eλnt]) را مختصرا E(t) بنامیم، نهایتا حل دستگاه غیر مزدوج را به این شکل داریم:x(t)=P E(t) P1 x(0)

مثالی از حل یک دستگاه معادلات دیفرانسیل مزدوج با استفاده از بردارهای ویژه

در مثال دستگاه معادلات دیفرانسیل مزدوج [x˙1x˙2]=[1302][x1x2] حل شده است.

حل یک دستگاه معادلات دیفرانسیل مزدوج با استفاده از بردارهای ویژه

مجددا با استفاده اکتاو اقدام به رسم نگارهٔ فاز ماتریس به همراه مسیر فاز چند نقطه می‌کنیم.

clear
[x1,x2]=meshgrid(-2:0.2:2,-2:0.2:2);
figure;
dx1= -1.*x1-3.*x2;
dx2= 0.*x1+2.*x2;
quiver(x1,x2,dx1,dx2);
hold on;

f=@(t,X) [-1*X(1)+-3*X(2);0*X(1)+2*X(2)];
a=[-2,-2,2,2];
b=[0.1,-0.25,0.2,-0.1];
for k=1:4
  [t,Sap]=ode45(f,[0 2],[a(k),b(k)]);
  plot(Sap(:,1),Sap(:,2),'LineWidth',2,'color','k')
end

v = linspace(-2,2,3)
w = -v
plot(v,w)
plot([-2,2],[0,0]);

axis equal;
xlabel('x1');
ylabel('x2');
title("Phase portrait for matrix A");
xlim([-2 2]);
ylim([-2 2]);

grid on;
Fig.1: نگارهٔ فازی که با ضرایب ماتریس A در مثال بالا ترسیم شده است.

توجه داشته باشید که اگر ماتریس A را پس از قطری شدن در نظر بگیریم (یعنی به ازای P1AP=[1002]) و نمودار متناظر با آن را رسم کنیم، چنین شکلی خواهیم داشت:

clear
[x1,x2]=meshgrid(-2:0.2:2,-2:0.2:2);
figure;
dx1= -1*x1-0.*x2;
dx2= 0.*x1+2*x2;
quiver(x1,x2,dx1,dx2);
hold on;

f=@(t,X) [-X(1);2*X(2)];
a=[0.2,-0.2,-0.2,0.2];
b=[2,-2,2,-2];
for k=1:4
  [t,Sap]=ode45(f,[0 2],[a(k),b(k)]);
  plot(Sap(:,1),Sap(:,2),'LineWidth',2,'color','k')
end

v = linspace(-2,2,3)
w = 0
plot([0,0],[2,-2])
plot([-2,2],[0,0]);

axis equal;
xlabel('y1');
ylabel('y2');
title("Phase portrait for matrix P-1AP");
xlim([-2 2]);
ylim([-2 2]);

grid on;
Fig.2: نگارهٔ فاز به ازای P1AP=[1002]

همانطور که مشاهده می‌شود از نگارهٔ ساده‌تر اخیر هم می‌توان به نگارهٔ پیچیده‌ی بالا رسید. به این صورت که اگر محورهای شکل Fig.2 را y در نظر بگیریم، تبدیل آن به محورهای x نمودار Fig.1، همانطور که در اثبات راهکار داشتیم، با استفاده از x=Py رخ می‌دهد. از آنجایی که نگارهٔ فاز Fig.2 از روی یک ماتریس قطری ساخته شده، می‌دانیم که بردارهای ویژه‌ی آن e1=[10] و e2=[01] هستند؛ روی شکل نمودار، مسیری که این دو می‌سازند با خطوط قرمز روی محورها مشخص شده است. حال بیاید تبدیل یافته‌ی بردارهای ویژه‌ی این نگاره را بیابیم. داریم: x=Py[1101][10]=[10] و x=Py[1101][01]=[11]همانطور که می‌بینید e1=[10] در نموداری که پایه‌هایش بردارهای ویژه نباشند به [10] بدل می‌شود و e2=[01] به [11]. پس بردارهای ویژه‌ی جدید را داریم. خطوط قرمز رنگ روی نگارهٔ Fig.1 هم‌راستا با همین بردارها ترسیم شده‌اند. سایر نقل و انتقالات از جمله تبدیل بردارها با استفاده از نگارهٔ ساده‌تر Fig.1 به Fig.2 نیز به راحتی با همین x=Py ممکن است. مثلا اگر یک نقطه نقطه‌ی تعادل باشد یا نقاط زینی این‌ها در حین تبدیل پایه‌ها تغییر نمی‌کنند.

الگو:تمرین توجه داریم که در اینجا می‌توانستیم به دو صورت بردار P را تشکیل دهیم، یکی با استفاده از P=[v1v2] و دیگری P=[v2v1]. ما در اینجا تعمدا P را به‌گونه‌ای ساختیم که جهت بردارها تغییر نکند زیرا همانطور که می‌دانیم P1x=y، و به ازای P1=[0111] ساخته شده از روی P=[1110]، چنین چیزی خواهیم داشت: P1x=[0111][x1x2]=[x2x1+x2]=[y1y2] یعنی x2 با محور y1 متناظر خواهد شد که مطلوب ما نیست.

شکل در صورت تعیین اشتباه محورها به همراه مسیرهای فاز (صرفا به جهت فلش‌ها توجه کنید که تغییر کرده است؛ اسم محورها و مسیرهای فاز را نادیده بگیرید)

الگو:پایان تمرین

جستارهای بیرونی

منابع

الگو:پانویس