رفتار سختی در ode45

هنگامی که برای ode45 خروجی در نظر نمی‌گیریم، در حالی که محاسبات اجرا می‌شود، همزمان نمودار نیز رسم می‌شود. نمودار از نقطه y=0.01 شروع می‌شودو با یک نرخ افزایشی ملایم رشد می‌کند تا اینکه مقدار T به 50 می‌رسد. در ان لحظه تا رسیدن به مقدار 1 رشد خیلی سریع اتفاق می‌افتد و نهاتا در مقدار 1 باقی می‌ماند.
مقدار δ را کوچکتر می کنیم تا میزان سختی این مسئله بیشتر احساس شود:

delta=0.0001;
ode45(F,[0 2/delta],delta,opts);

خروجی مانند شکل قبل است با این تفاوت که زمان خیلی بیشتری در این مرحله نیاز است. اگر خسته شدید می‌توانید از دکمه stop در سمت چپ تصویر استفاده کنید.
این مسئله در ابتدا سخت نیست، ولی به محض اینکه به حالت پایدار می‌رسد سخت می‌شود. هر جوابی نزدیک به y(t)=1 به سرعت به جواب مسئله نزدیک می‌شود.
برمی‌گردیم به مثال! این بار با یکی از odeهای MATLAB که نام آن با “s” برای "stiff" ختم می‌شود استفاده می‌کنیم.

delta=0.0001;
ode23s(F,[0 2/delta],delta,opts);

می‌توانید ببینید که ode23s گام‌های کمتری برای حل معادله دیفرانسیل نسبت به ode45 می‌گیرد. در واقع این مثال یک مسئله ساده برای یک حل‌کننده "سخت" است. سختی حتی بر روی خروجی گرافیکی نیز تاثیر می‌گذارد. در صورتی که بخواهید خروجی‌ها را چاپ کنید، کیفیت ode23s خیلی بهتر از ode45 خواهد شد.
یک نکته جالب که در مورد مسئله شعله آتش وجود دارد، وجود تابع W لمبارت، W(z)، است. معادله مذکور یک معادله جداشدنی است. تابع y بر حسب t بصورت زیر بیان می‌شود:

این معادله می‌تواند بر حسب y نوشته شود. فرم دقیق تحلیلی مدل شعله بصورت زیر خواهد بود:

که در این معادله a=1/delta-1 . تابع W(z)، تابع W لمبارت است که در واقع حل معادله زیر می‌باشد:

با استفاده از MATLAB و جعبه ابزار Symbolic متصل به Maple، معادله دیفرانسیل مربوط به مدل شعله را حل می‌کنیم تا به یک جواب دقیق برسیم و آن را رسم می‌کنیم.

y = dsolve('Dy = y^2 - y^3','y(0) = 1/100');
y = simplify(y);
pretty(y)
ezplot(y,0,200)

1
-------------------------------
lambertw(0, exp(99 - t) 99) + 1

References
Robert M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, and D. E. Knuth, “ On the Lambert W Function,” Advances in Computational Mathematics, Volume 5, 1996, pp. 329-359.

نویسنده: (#جبار_کمالی )
#For_intermediate , #ODE, #Stiffness
@MatlabTips

🔵کنترل تولید اعداد تصادفی🔵
سطح پیچیدگی:🌕🌑🌑🌑🌑
پیش‌نیاز:کدنویسی پایه
برای:همه
📝 هدف از این پست: (استفاده حرفه ای تر از تابع رندم و یک گذار به کاربردهای آن در رمزنگاری)


ممکن است برای شما این سوال پیش آمده باشد که چرا ما اینقدر بر روی تولید اعداد تصادفی تاکید میکنیم و اصلا چه لزومی برای این کار وجود دارد. ساده ترین جواب این است که بسیاری از دوستان از اعداد تصادفی به طرق اشتباهی در برنامه هایشان استفاده می کنند و از آنجایی که روند ایجاد اعداد تصادفی چیزی نیست که به راحتی بتوان دیباگ کرد(به خاطر اینکه ما انتظار رفتار خاصی از یک دنباله تصادفی نداریم) یافتن منبع اشتباه در برنامه بسیار دشوار می شود. برای حل این مشکل متلب تابعی به نام rng قرار داده است که نوع تولید اعداد تصادفی را مشخص کرده و هر بار دانه اولیه را هم میگیرد. به این ترتیب می توانید الگوریتم خودتان را با یک دنباله تکرار شدنی از اعداد تصادفی تست کنید و سپس قسمتهایی از کد را که فکر می کنید مشکل دارند تغییر دهید. برای شروع در خط فرمان تابع زیر را وارد کنید تا مقادیر را ببینید

>> t = rng      %Stores current settings of rng in t
t = 

     Type: 'twister'
     Seed: 0
    State: [625x1 uint32]

همانطور که می بینید نوع ساخت عدد تصادفی Mersenne Twister است و دانه آن برابر با صفر است و وضعیت آن هم در یک ماتریس ذخیره شده است. حالا دستور rand را وارد میکنیم تا نتیجه را ببینیم:

>> rand(1,4)     %Generates a 1x4 random matrix with default initial
                 %state
ans =

    0.8147    0.9058    0.1270    0.9134
>> rand(1,4)     %Generates a 1x4 random matrix
                 %(different from initial random matrix)
ans =

    0.6324    0.0975    0.2785    0.5469

سپس تولید اعداد شبه تصادفی را به حالت اول بر میگردانیم و دوباره تابع rand را اجرا میکنیم

>> rng(t)        %Restore settings of rng
>> rand(1,4)     %New random matrix same as initial random matrix
ans =

    0.8147    0.9058    0.1270    0.9134

همانطور که می بینید دوباره همان اعداد تولید شدند! این برای دیباگ کردن برنامه هایی که از اعداد رندم استفاده میکنند عالی است!
نسخه کنونی متلب از الگوریتم Mersenne Twister برای تولید اعداد شبه تصادفی استفاده میکند[1]. الگوریتم های دیگری هم برای این کار وجود دارند که می توانید در help متلب با تایپ کردن rng بیابید. (مانند SIMD fast Mersenne Twister، Combined Multiple Recursive و Multiplicative Lagged Fibonacci)
هر بار که متلب ریستارت می شود. دانه آن به مقدار قبلی (یعنی صفر) ست می شود!( اگر باور نمی کنید متلب را باز کنید تابع Rand را اجرا کنید مقدار را به خاطر بسپارید و سپس متلب را ببندید و دوباره باز کنید و باز هم Rand را بزنید) برای اینکه این مشکل را حل کنید دستور زیر را وارد کنید:

>> rng('shuffle')     %Seed rng on basis of current time
>> rng
ans = 

     Type: 'twister'
     Seed: 385954034
    State: [625x1 uint32]
>> rand(1,4)          %New random matrix produced with new seed
ans =

    0.5706    0.9136    0.4043    0.8423

این دستور باعث می شد دانه الگوریتم به زمان تنظیم آن ست شود(یک راه حل ساده و هوشمندانه!).
در قسمت قبل شما را تشویق کردم که در مورد تولید اعداد واقعا تصادفی فکر کنید! از آنجا که کامپیوتر های کنونی ماشین های محاسبه متعین(deterministic ) ایی هستند نمی توان هیچ الگوریتمی برای تولید اعداد تصادفی ایجاد کرد. مهمترین روش تولید اعداد تصادفی استفاده از یک فرآیند فیزیکی مانند نویز حرارتی، اثر فوتوالکتریک یا پدیده ای کوانتومی است. این فرآیند ها غیر قابل پیش بینی هستند[2] و بنابراین پایه ای برای تولید اعداد واقعا تصادفی هستند. مثلا میزان گرمای حرارتی پردازنده یا حالت کوانتومی بخشی از حافظه می تواند مورد استفاده قرار گیرد.
سیستم عامل مسیری برای تولید این اعداد را فراهم میسازد. تولید این اعداد از ویندوز 2000 دچار مشکلاتی بود که در سال 2008 برطرف شد[3]. در سیستم عامل های یونکیسی هم از /dev/random استفاده می شود که بسیار مقاوم تر(robust) از نسخه های ویندوزی است. اعداد تصادفی واقعی بیشترین کاربرد را در رمزنگاری و شبیه سازی (مانند Monte carlo) دارند. از آنجا که حمله های مبتنی بر تولید اعداد تصادفی بسیار مهلک هستند تولید اعداد "واقعا" تصادفی بسیار مهم است. سایتی به نام random.org ادعا می کند که اعداد واقعا تصادفی بر اساس نویز اتمسفری تولید میکند که برای استفاده از آن در متلب یک واسط وب طراحی شده است که من تغییراتی در آن داده ام تا قابل استفاده باشد. در زیر می توانید دانلود کنید

استفاده از این اسکریپت بسیار ساده است. چند مورد را ببینید: ( می توانید آن را به path اضافه کنید تا همیشه داشته باشیدش)

>>truerand %create a true random integer from 1 to 100

ans =

    60
>> truerand(3,6) % create a matrix of 3x6 of true random numbers

ans =

    77    67    30    25    94    86
    51    25    76    86     4    20
    26    95    75    57    22    48

به خاطر داشته باشید که هر بار که این تابع را صدا می زنید یک HTTPS GET انجام می شود که البته زمانبر است. پس در صورتی که خارج از ایران هستید یا نیاز شدید دارید می توانید از آن استفاده کنید.
—-------------------------------------------------------------
[1] https://en.wikipedia.org/wiki/Mersenne_Twister
[2] Pironio, Stefano, Antonio Acín, Serge Massar, A. Boyer de La Giroday, Dzimitry N. Matsukevich, Peter Maunz, Steven Olmschenk et al. "Random numbers certified by Bell’s theorem." Nature 464, no. 7291 (2010): 1021-1024.
[3] https://www.computerworld.com/article/2539986/security0/microsoft-confirms-that-xp-contains-random-number-generator-bug.html


@MatlabTips
#For_beginner
#Rand #Truerand #Cryptography
نویسنده:(#A_2)

🔵درباره dsolve و ode🔵
سطح پیچیدگی:🌓🌑🌑🌑🌑
پیش‌نیاز: (ندارد)
برای:(همه)
📝هدف از این پست: (این پست کاربرد و تفاوت دو دستور فوق العاده کاربردی در حوزه معادلات دیفرانسیل در متلب بیان میکند)
🔵dsolve🔵
دستور dsolve مربوط به کتابخانه سمبولیک ریاضی متلب است. پس نتیجتا خروجی هم بصورت سمبولیک خواهد بود. سمبولیک بودن چیزی جدای نمایشی مانند float-point یا نمایش های عددی دیگر است. تو حالت سمبولیک بیشتر با نماد سر و کار داریم تا عدد. شما می توانید یک معادله یا دستگاهی از معادلات مختلف را به راحتی با دستور dsolve حل کنید. خروجی به یکی از سه صورت زیر خواهد بود:
1⃣خروجی به صورت صریح بر حسب توابع مقدماتی(مثل چندجمله ای ها، مثلثاتی و ...)تولید می شود.
2⃣ خروجی به صورت صریح بر حسب توابع خاص(توابعی که خود جواب یک معادله دیفرنسیل معروف هستند مثل gamma) است.
3⃣ خروجی به صورت ضمنی تولید می شود.
بزرگترین مزیت دستور dsolve این است که یک فرمول بر میگرداند(البته نه همیشه) و البته خیلی از معادلات هم جوری حل میکند که هیچکس را خوش نمی اید.
🔵ode🔵
این نوع دستور دیگر مانندdsolve از کتابخانه سمبولیک نیست و یک ام فایله. این توابع معادله را با استفاده از الگوریتم هایی که در تکست های دانشگاهی (مانند رانگ-کوتا) به صورت عددی حل میکنند. برای مثال الگوریتم ode23 بر اساس مقالات Bogacki و Shampine است. عدد 23 در این تابع به دلیل وجود دو الگوریتم single-step(تشریح این الگوریتم در اسکوپ این پست نیست) است که یکی از آنها درجه 2 و دیگری درجه 3 است. مشابه این تابع، ode23s هست که در آن s نماینده کلمه "stiffness" است که پستی مخصوص به همین لغت قبلا گذاشته شد.
نویسنده:(#جبار_کمالی)
#For_all , #ODE
@MatlabTips

🔵نمونه گیری تصادفی🔵
سطح پیچیدگی:🌕🌕🌑🌑🌑
پیش‌نیاز:کدنویسی پایه
برای:همه
📝 هدف از این پست: (استفاده از امکانات جدید متلب برای نمونه گیری تصادفی و کاربرد آن در حالا مختلف)


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

out = randsample([3,5,1,6,7], 3)

در این قطعه کوچک برنامه شما سه عدد به صورت تصادفی از لیست اعداد [3,5,1,6,7] انتخاب می کنید. حالا ممکن است شرایط قدری پیچیده باشد یعنی شما می خواهید انتخاب های شما با جایگذاری باشند. مانند یک کیسه توپ که پس از انتخاب دوباره آن را به کیسه برگردانید تا شانس دوباره ای به آن داده باشید. در این حالت کافی است آرگومان سوم(replacement ) را true کنید.

out = randsample([3,5,1,6,7], 3, 'true')

حالت های پیچیده تر هم ممکن است پیش بیاید! من در این جا یک مثال کاربردی ذکر می کنم: در الگوریتم کلونی مورچگان در قسمتی از الگوریتم(مهمترین قسمت آن) مورچه باید از بین چندین مسیر از نودی که در آن قرار دارد یکی را به صورت تصادفی انتخاب کند. . یک نما از این حالت را در شکل زیر می بینید:

اما مساله این است که احتمال انتخاب مسیر ها با هم متفاوت است و بنابراین این احتمال را باید در انتخابمان لحاظ کنیم. فرض کنید پس از قرار دادن در فرمول (به الگوریتم کلونی مورچه مراجعه کنید) احتمال مسیر ها به صورت یک لیست 5 تایی به صورت [0.2, 0.1, 0.1, 0.05, 0.3, 0.05, 0.12, 0.08] به دست آمدند که به ترتیب متناظر با حرکت به سمت نود های [3,5,1,6,7] هستند. حالا چون ما می خواهیم فقط یک مسیر را انتخاب کنیم آرگومان چهارم که احتمال مسیر هاست را به صورت زیر به خورد تابع می دهیم:(وقتی یکی انتخاب می کنیم قاعدتا با جایگذاری یا بدون آن مهم نیست)
‍‍‍‍‍‍‍

SelecteNode = randsample([3,5,1,6,7], 1, “true”,[0.2, 0.1, 0.1, 0.05, 0.3, 0.05, 0.12, 0.08] )

فراموش نکنید که لزومی ندارد حتما یک زیر مجموعه از نمونه ها را انتخاب کنید! شما می توانید بر اساس توزیع احتمالی که دارید نمونه تولید کنید. مثلا فرض کنید که می خواهیم یک دنباله 48 تایی از دنباله DNA را تولید کنیم با این دانش که از چهار پروتئین سازنده DNA یعنی (Adenine, Thymine, Guanine , Cytosine) احتمال حضور آن هریک در زنجیره ی مارپیچی به صورت زیر است:

R = randsample('ACGT',48,'true',[0.15 0.35 0.35 0.15])

با همین تابع شما می توانید مقدار زیادی کد، وقت و اعصاب از خودتان را ذخیره کنید!
در قسمت بعدی نشان می دهیم که چگونه می توان همین کار را برای یک تابع احتمال پیوسته انجام داد.

@MatlabTips
#For_All
#Ransample #Rand #Ant_colony #Random_sample
نویسنده:(#A_2)

🔵نگاهی کیفی به معادلات دیفرانسیل🔵
سطح پیچیدگی:🌕🌓🌑🌑🌑
پیش‌نیاز: برای همه
📝هدف از این پست:( آشنایی با دستور quiver و رفتار میدان های برداری)

می توان جواب یک معادله دیفرانسیل را بدون اینکه آن را حل کنیم بصورت جئومتریکی بررسی کنیم. یکی از این روش ها استفاده از میدان های برداری است.
میدان برداری با رسم خط کوتاهی که از نقطه (t,y) عبور میکند و شیب آن f(t,y) است. جواب ها، یا منحنی‌های انتگرالی، این ویژگی را دارند که در هر نقطه مماس بر میدان‌ برداری است و بنابراین یک ارزیابی کلی و کیفی از حل یک معادله دیفرانسیل از میدان برداری استنتاج می شود.
برای معادلات ساده، می توانید میدان برداری را دستی رسم کنید، اما Matlab برای هر معادله درجه اولی قادرست میدان برداری را رسم کند.
در زیر یک مثال برای معادله دیفرانسیل خطی آورده شده است.

دستوری که متلب برای رسم میدان برداری استفاده میکند quiver است. در کنار این دستور همواره meshgrid وجود دارد.

t = -2:0.1:3;
y = -1:0.1:2;
[T,Y]= meshgrid(t,y);
S = Y.^2 + T;
L = sqrt(1 + S.^2);
quiver(T,Y,1./L,S./L,0.5)
axis tight
title 'superimpose of direction field and y''(t) = y^2 + t'
axis([-2 3 -1 2])

کاری که این کد انجام میدهد مجموعه ای نقاط (گرید) میسازد و این نقاط را در معادله دیفرانسیل (شیب هر نقطه) قرار میدهد. بردارهای (1,s) شیب مطلوب را به ما میدهد ولی با طول های متفاوت. به همین منظور نرمال کرده ایم.