برای دیدن فهرست باقی نوشته‌های این مجموعه به قسمت اول (کلیات) بروید
در تحلیل اجزا محدود ما با ماتریس‌ها خیلی سر و کار داریم. هر گره به مساله‌ی ما تعدادی درجه‌ی آزادی اضافه میکنه و در نهایت ما نیاز داریم تا با ماتریسهای خیلی بزرگی کار کنیم. بعنوان مثال این مدل المان محدود رو در نظر بگیرید:
این مدل رو به‌صورت یک جدول (Grid) سه بعدی در نظر میگیریم که ابعادش N هست و در هر بعد تعداد N عدد گره داره. برای این مثال، N برابر 9 هست. در این مدل تعداد کل گره ها برابر N به توان 3 میشه که برای این مثال میشه 9بتوان3 برابر 729 تا گره. و اگر برای هر گره سه درجه آزادی غیر چرخشی در نظر بگیریم، با ضرب 729 در 3 برابر حدود 2200تا درجه آزادی داریم. تعداد المانهامون میشه N-1بتوان 3 که N=9 میشه 512 عنصر خواهیم داشت که هر کدوم 8 تا گره رو به هم متصل میکنند و المانهای مکعبی هستند. و هر گره هم سه درجه ازادی پس یعنی ماتریس سختی هر عنصر ماتریس مربعی خواهد بود که ابعادش 24 در 24هست. یعنی حداکثر میتونه برابر 24*24 برابر 576تا عضو غیر صفر داشته باشه. این ماتریس ها رو باید در ماتریس سختی کل سوار کنیم. با ضرب 576 در تعداد عناصرمون که 512 هست، پس حداکثرِ ممکن حدود 295K عضو غیر صفر داریم که در ماتریس سختی اصلی سوار خواهند شد. البته این مقدار دقیق نیست یجورایی تئوری هست و مقدار اصلی به نظرم از این کمتر میشه. میخوام بگم این یعنی حداکثر فقط 6 درصد اعضای ماتریس سختی کل غیر صفر هستن و باقیشون صفر هستن. هرچی مدل بزرگتر بشه این درصد کمتر میشه، مثلا اگر بجای 9، ما مقدار N رو برابر 70 در نظر میگرفتیم که عدد دور از ذهنی هم نیست این مقادیر این طور تغییر میکنند:
در ادامه مثال بالا رو با فرض N=70 دنبال میکنیم. در این حالت اگر ماتریس سختی کل رو با استفاده از آرایه‌ی (Array) دوبعدی در زبان سی‌شارپ معرفی کنیم، خواهیم داشت:
double[,] matrix = new double[1000000, 1000000]
با این کد یک نمونه از شی آرایه ساخته خواهد شد و رانتایم سعی خواهد کرد تا یک فضایی رو برای این آرایه در نظر بگیره. نوع داده‌ای که ما استفاده کردیم double هست. در خیلی از زبان‌های برنامه‌نویسی این نوع داده وجود داره، و یک تعریف کلی و استاندارد داره که از اینجا اطلاعات بیشتر در دسترس هست: لینک https://en.wikipedia.org/wiki/Double-precision_floating-point_format
طبق تعریف، هر نمونه از این نوع داده ۸ بایت فضا نیاز داره و بطور کلی ما ماتریسی با ابعاد ۱م. (بخوانید ۱ میلیون) در ۱م. داریم. تعداد اعضای این ماتری برابر ۱م. به‌توان ۲ هست یعنی ۱۰۰۰میلیارد عضو. یا ۱۰۰۰ب. (بخوانید ۱ میلیارد، یا بیلیون به قول خارجیا). هر عضو ۸ بایت فضا نیاز داره و با ضرب ۸ بایت در ۱۰۰۰میلیارد متوجه میشیم که ۸ ترابایت فضا نیاز خواهیم داشت، اون هم در RAM و نه هارد. وقتی این دستور رو بخواهیم با زبان سی شارپ اجرا کنیم یک خطا دریافت میکنیم به عنوان اینکه رم‌نداریم آقا چه خبره؟! حداقل با کامپیوترهای امروزی (الان سال 1400 شمسی!) که اغلب کامپیوترهای شخصی بین 4 الی 32 گیگ RAM دارند این روش قابل پیاده سازی نیست. بعلاوه برای یک ضرب ساده این ماتریس در یک بردار یه تعداد حدودا 2T (دو تریلیون یا دو ترا) ضرب و جمع باید انجام بشه که برای کامپیوترها میتونه زمانبر باشه. حل دستگاهش که دیگه بماند، چون خیلی خیلی بیشتر از یک ضرب ساده عملیات خواهد داشت.
به‌خاطر همین محدودیت، روش دیگه ای ابداع شد که ماتریس خَلوَت، تُنُک یا Sparse نام داره. (دقیقا عین اعصاب من که مثل موهای یکی از بهترین دوستام خلوته ). در این روش ما میاییم بجای اینکه یک ماتریس خیلی بزرگ در نظر بگیریم که تقریبا ۹۹ درصد اعضاش صفر هستند، فقط اعضای غیر صفر رو در نظر میگیریم. به‌عنوان‌مثال یک لیست که شامل سطر و ستون و مقدا اعضای غیر صفر هست. در این حالت حافظه مورد نیاز خیلی کمتر میشه، یعنی برای هر عضو غیر صفر دو متغیر نوع صحیح (integer) برای ذخیره شماره سطر و ستون و یک متغیر نوع اعشاری (double یا float64) برای مقدار خواهیم داشت که میشه جمعا 4بعلاوه4بعلاوه8 بایت برابر 16 بایت. پس بجای 8 ترابایت رم، به مقدار 190M (تعداد اعضای غیر صفر یا NoNZero یا nnz) ضربدر 16بایت، برابر حدودا 3 گیگابایت حافظه رم نیاز داریم که خوشبختانه روی اغلب کامپیوترها وجود داره.
برای همین مثال ابعاد بردارهای نیرو و تغییر مکان برابر 1M در 1 خواهد بود که هر کدام به 8مگابایت حافظه برای ذخیره شدن نیاز دارند. چون این مقدار نسبتا کم هست، و ما دغدغه‌ی فراهم کردن ۸ مگابایت رم رو نداریم. پس دیگه از روش ماتریس خلوت برای بردارها استفاده نمی‌کنیم چون استفاده از ماتریس خلوت هزینه داره و هزینه‌اش پیچیده شدن برنامه هست. پس برای بردارها هم عناصر صفر و هم غیر صفر را ذخیره میکنیم. این روش رو اصطلاحا ماتریس متراکم یا Dense Matrix می‌نامیم. ضمنا طبق تعریف به ماتریسی که تعداد ستون برابر 1 و تعداد سطر متغیر دارد بردار میگوییم. پس بردار به طول 1M منظورمان ماتریسی با ابعاد یک ستون در 1M سطر هست.
روشهای مختلفی برای ذخیره سازی ماتریس خلوت وجود داره. یکی از ساده‌ترین‌هاش روش "فهرست مختصات" (Coordination list) یا به اختصار COO، همینی بود که بالا در موردش صحبت کردیم. روش دیگه "ردیف خلوت فشرده" (Compressed Sparse Row) یا به اختصار CSR یا CRS یا YALE نام داره که کمی پیچیدگی بیشتری داره. فهرستی از روشها اینجا قابل مشاهده هست.
https://en.wikipedia.org/wiki/Sparse_matrix
در روش CSR مانند روش COO از سه‌تایی مرتب استفاده میکنیم با این تفاوت که بجای row تعداد عناصر غیر صفر در اون سطر رو مینویسیم. عکس زیر بهتر توضیح میده:
این رو خوب توضیح ندادن، خودتون سرچ منید یا برای اطلاعات جامع تر لطفا به این پیوندها مراجعه کنید
https://blog.faradars.org/%D9%85%D8%A7%D8%AA%D8%B1%DB%8C%D8%B3-%D8%AE%D9%84%D9%88%D8%AA-%D8%AF%D8%B1-%D8%B1%DB%8C%D8%A7%D8%B6%DB%8C%D8%A7%D8%AA-%D9%88-%D8%B3%D8%A7%D8%AE%D8%AA%D9%85%D8%A7%D9%86-%D8%AF%D8%A7%D8%AF%D9%87/
https://developpaper.com/instructions-for-using-python-scipy-sparse-matrix/
با استفاده از ماتریس خَلوَت یا Sparse میشه ماتریسِ سختیِ مدل‌هایِ بزرگ رو در حافظه ذخیره کرد. بعنوان مثال توانستیم ماتریس سختی یک مدل با 1M درجه آزادی رو بجای 8ترابایت فضا، در 3گیگابایت فضا ذخیره کنیم. برای ذخیره ماتریسهای بزرگ مانند ماتریس سختی، جرم و میرایی از ماتریس خلوت استفاده میکنیم. ولی بردار نیرو و تغییر مکان رو بصورت ماتریس متراکم یا (Dense) در نظر میگیریم.
عملیات ضرب و جمع بین ماتریسهای خلوت به سادگی ماتریس متراکم نیست و پیچیده تره. بهتره مسایل مربوط به ماتریسهای خلوت بصورت جدا از مساله اصلی توسعه پیدا کنه. یعنی توسط یک گروه دیگه.
کتابخانه ی BFE از کتابخانه ای دیگر بنام CSParse.NET استفاده میکنه که با زبان C# نوسته شده. خود این تابخانه در اصل برگردانی هست از برنامه قدیمیتر دیگری بنام CSparse که آقایی بنام Tim Davis اون رو نوشتن به زبان C. این یک کتابخانه نسبتا کامل و بسته آماده هست که مسایل مربوط به ماتریس خلوت، اعم از ساخت، ترانهاده کردن، ضرب، جمع، تجزیه و حل دستگاه معادلات خلوت رو انجام میده و لازم نبود که ما این کدها رو خودمون بنویسیم. در قسمت حل معادلات که در ادامه میاد، کمی با جزییات محاسبات ماتریس خلوت آشنا خواهید شد، بعلاوه اینکه اگر این کدها رو خودمون مینوشتیم بعدها برای نگهداری (Maintenance) و باگزدایی (Debugging) مشکلات عدیده ای پیدا میشد چون آدم یادش میره چی نوشته... به این ترتیب با استفاده از github کار گروهی میکنیم.
ضمنا در مورد ماتریس Sparse قسمت اعظم اطلاعاتم از کتابی به‌نام Sparse Matrix (لینک) نوشته‌ی آقای سرخیو پیزانسکی به دست اومده. به شما هم پیشنهاد میکنم اگر تمایل و نیاز داشتید کتابش رو نگاهی بی‌اندازید.
فرض کنید ما یک دستگاه معدلات داریم در فرم ماتریسی:
F=K.U
در این رابطه F و U به ترتیب بردار تغییر مکان و نیرو هستند (تعداد ستون بردارها همیشه 1 هست و تعداد سطر متفاوت دارند) و ماتریس K همون ماتریس سختی مون هست. دقت داشته باشید ماتریس K به خودی خود تکینه (Singular) هست و این به این معنیه که از نظر جبری معکوس نداره و دترمینانش برابر 0 هست. اما بعد از اینکه سطرها و ستونهای مربوط به درجات آزادی رو ازش حذف کنیم مشکل تکینه بودن رفع میشه و دستگاه قابل حل میشه. خوب برگردیم سر اصل مطلب، همون دستگاه معادلاتمون. فرضا ابعاد ماتریس K برابر تعداد سطرهای F و U و برابر n هستند (تعداد درجات ازادی غیر مقید). توی این حالت ما ماتریس K رو داریم و برایمان مجهول نیست، ماتریس F رو هم داریم و برایمان مجهول نیست، مقادیرش نیروهای وارده بر درجات ازادی غیر مقید هست. بردار U در این دستگاه مجهول هست. خوب اولین چیزی که به ذهنمون میرسه اینه که K رو معکوس میکنیم و در F ضرب میکنیم و در نهایت ماتریس U را میبایم . بیایید ببینیم اگر بخواهیم این روش رو پیاده کنیم باید چکار کنیم. فرض کنید تابعی داریم بنام Invert که یک ماتریس رو بعنوان ورودی میگیره و معکوس اون رو بعنوان خروجی بهمون میده. و یک تابع بنام Multiply برای ضرب ماتریس در بردار
var invK = Invert(K); var U = Multiply( invK,F);
با این روش میتونیم U رو بدست بیاریم، منتهای مراتب مشکلاتی وجود داره. مهمترین این مشکل ها اینه که معکوس یک ماتریس خلوت الزاما یک ماتریس خلوت نیست. توی قسمت قبل دیدیم چطور ماتریس K دارای تعداد زیادی عنصر صفر هست و در یک مثال هم دیدیم نسبت عناصر غیر صفر به کل عناصر در ماتریسمون حدود دو صدم درصد میشه و میتونیم بجای 8 ترابایت حافظه در 3 گیگابایت حافظه ذخیره اش کنیم. حالا اگر این ماتریس رو معکوس کنیم و در ماتریس معکوس شده نسبت اعضای غیر صفرش مثلا 50 برابر بشه، یعنی بشه 1 درصد، اون وقت به مشکل میخوریم چون 50 برابر بیشتر حافظه RAM نیاز داریم. اون هم فقط برای ذخیره کردن ماتریس توی حافظه. انجام اعمال ضرب و تقسیم روی این تعداد عدد بماند. پیشنهاد میکنم این مقاله رو بخونید: https://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/ که در مورد حل دستگاه معادلات Sparse هست. روشهایی وجود داره که میشه بدون معکوس کردن ماتریس K دستگاه رو حل کرد. این روشها به دو دسته تقسیم میشن: روشهای تکراری (Itterative) و روشهای مستقیم (Direct).
ساده ترین این روشها، روش حذفی گاوس هست که با استفاده از جمع و تفریقِ سطرها، اعضای غیر صفر ماترسی K رو حذف میکنیم.https://fa.wikipedia.org/wiki/%D8%AD%D8%B0%D9%81_%DA%AF%D8%A7%D9%88%D8%B3%DB%8C که این روش زیاد بطور مستقیم استفاده نمیشه.
یک دسته روشهایی هستند که ماتریس K رو به ضرب دو یا سه ماتریس در هم تجزیه میکنند. مثل روش LU یا LDL یا Cholesky که در زمره ی روشهای مستقیم قرار میگیرن. با این روشها هم الزاما ماتریس بعد از تجزیه بصورت Sparse باقی نمیمونه، ولی با یک سری ملاحظات میشه تعداد اعضای غیر صفر رو کم کرد.
روشهایی مثل گاوس-سایدل یا گرادیان مزدوج (Conjugate Gradient) در زمره ی روشهای تکراری قرار میگیرن. یک روش مانند روش حذفی گاوس در نگاه اول چندان پیچیده نیست. ابتدا باید یک سطر رو در نظر بگیریم و اصطلاحا اون رو پاشنه (Pivot) قرار بدیم و ضریبی از اون رو به سطرهای دیگه اضافه کنیم به نحوی که تمام ستون صفر بشه به غیر از پاشنه. این روش برای ماتریسهای متراکم (Dense) خیلی راحت قابل پیاده سازی هست. ولی برای ماتریس های خلوت به این راحتی ها نیست. با هر بار عملیات سطری مقدماتی الگوی غیر صفر ماتریس ما تغییر میکنه و در نتیجه حافظه ی مورد نیاز هم تغییر میکنه. برای این کار اول میان تحلیل سمبلیک (Symbolic Analysis) انجام میدن به این معنی که بدون در نظر گرفتن مقادیر، ابتدا ساختار غیر صفر ماتریس نهایی رو محاسبه میکنن، و حافظه رو تخصیص میدن بعد مقادیر رو توش پر میکنن و عملیات سطری مقدماتی رو انجام میدن و در نهایت جواب رو بدست میارن. به مقدار فضای اضافه شده نسبت به الگوی غیر صفر اولیه اصطلاحا Fill In میگن. در صورتی که ماتریس رو بدون ترتیب بندی (reordering) سطر و ستونها تحلیل سیمبولیک کنیم Fill In خیلی بالا میره. ولی در صورتی که دوباره ترتیب بندی کینم Fill In به مقدار کمینه ی خودش نزدیک تر میشه. دقت کنید یافتن بهترین ترتیب که مقدار fill رو به حداقل برسونه یک مساله ی از اون هاست که نمیشه راحت حلش کرد ولی با روشهایی که دم دست داریم میشه بهش حداقل مقداری نزدیک شد. این روشها پیچیده هستند و رابته تنگاتنگی با علم تحلیل گراف دارن و از الگوریتمهاش استفاده میکنن.
طبق مستندات روش چولسکای هم مشابه همین روش حذفی گاوس هست یعنی بطر مشابه تحلیل سیمبولیک میشه.
خوشبختانه اون موقع که ما شروع به کار کردیم کتابخانه ای بنام CSparse.NET موجود بود که تمتا این جزییات رو پیاده کرده و یک راه حل حاضر و اماده در اختیار ما قرار داد. تغریبا تنها وابستگی ای که کتابخانه BFE داره به همین کتابخانه CSparse.NET هست. این کتابخونه با اینکه با زبان C# نوشته شده سرعت مناسبی هم داره و طبق آزمایشا هایی که نویسنده اش انجام داد سرعتی نزدیک به کتابخانه اصلی Csparse داره (در برخی موارد تا 80 درصد سرعت اون رو داره و اون با زبان C نوشته شده). الگوریتمهاش هم خیلی خوبه و نسبتا کامله. مهمترین نکته اش اینه که پشتیبانی داره و میشه از نویسنده اصلی کمک گرفت. این خیلی مهمه...
