طرق Runge-Kutta المباشرة. قائمة المصادر المستخدمة

07.05.2019

دعونا نعطي معادلة تفاضلية من الدرجة الأولى

مع الحالة الأولية

ص(س 0) = ص 0.

دعنا نختار الخطوة h ونقدم الترميز التالي:

س ط = س 0 + ط . ح و ص أنا = ذ(س أنا) ,

حيث أنا = 0، 1، 2، ...

س ط - العقد الشبكة،

y i - قيمة الدالة المتكاملة في العقد .

يتم حل المعادلة التفاضلية بشكل مشابه للطرق الموضحة أعلاه. الفرق هو أن الخطوة مقسمة إلى 4 أجزاء.

وفقا لطريقة Runge-Kutta من الدرجة الرابعة، القيم المتعاقبة ذ أنا الوظيفة المطلوبة ذيتم تحديدها بواسطة الصيغة:

, ط = 0، 1، 2، ...

والأرقام ك 1 (ط) , ك 2 (ط) , ك 3 (ط) , ك 4 (ط)في كل خطوة يتم حسابها باستخدام الصيغ:

هذه طريقة واضحة من أربع خطوات بدقة من الدرجة الرابعة.

تعتبر طرق Runge–Kutta سهلة البرمجة وتتمتع بقدر كبير من الدقة والثبات في التعامل مع نطاق واسع من المشكلات.

ويبين الشكل 6 مخططا انسيابيا لهذا الإجراء رونج (X0، XK، Y0، N، Y)لحل مسألة كوشي باستخدام طريقة رونج-كوتا الموضحة أعلاه.


ط = 0، …، ن-1

K1 = ح * و(س، يي)

K2 = ح * F(س + ح/2، يي + K1 / 2)

K3 = ح * و(س + ح/2، يي + K2 / 2)

K4 = ح * F(س + ح، يي + K3)

ك = (ك1 + 2*ك2 + 2*ك3 + ك4) / 6

الشكل 6 - مخطط انسيابي لإجراء RUNGE

يوضح الشكل 7 رسمًا تخطيطيًا لخوارزمية البرنامج الرئيسية لحل مشكلة كوشي والحصول على نتائج بعدد ثابت من مقاطع القسم N. في البرنامج الرئيسي، يسمى الإجراء رونج (X0، XK، Y0، N، Y)،حساب قيمة الوظيفة المطلوبة ذ يفي النقاط س يبطريقة رونج-كوتا.

البيانات الأوليةفي هذه المشكلة هي:

X0، XK - القيم الأولية والنهائية للمتغير المستقل؛

Y0 - القيمة ص 0من الحالة الأولية ص(س 0) = ص 0;

N - عدد أجزاء القسم.

نتائج البرنامجيتم عرضها في عمودين:

X - مجموعة من قيم عقدة الشبكة؛

Y - مجموعة قيم الحل المطلوب في عقد الشبكة المقابلة.


الإدخال X0، XK، Y0، N


رونج (X0،XK،Y0،N،Y)


ط=0…ن


الإخراج X، Y ط


الشكل 7 - رسم تخطيطي لخوارزمية البرنامج الرئيسية لحل مشكلة كوشي مع عدد ثابت من مقاطع القسم N


حل المعادلات التفاضلية في MathCad

الشكل 8 - مثال على حل معادلة تفاضلية باستخدام الطريقة

أوامر Runge-Kutta 4 في بيئة MathCad.

وظائف الرسوم البيانية في بيئة Visual Basic

تعد مهمة جدولة الدوال وإنشاء الرسوم البيانية الخاصة بها إحدى المهام الرئيسية في عملية حل المعادلات التفاضلية العادية. دعونا نفكر في هذه المهمة بمزيد من التفصيل.

مهمة.

أنشئ رسمًا بيانيًا للدالة y=sin(x) على القطعة. خذ خطوة الجدولة التي تساوي h.

حل.

لرسم رسم بياني للدالة في بيئة Visual Basic، من المناسب استخدام بعض المكونات الرسومية.

الشكل 9 - موقع المكونات الرئيسية في النافذة العامة



يتم استخدام مكون Picture Box() كحاوية للتخطيط. وهي عبارة عن مصفوفة من النقاط (البكسلات)، ويمكن التحكم في لون كل نقطة على حدة. يتم تحديد إحداثيات أي نقطة بواسطة زوج من الأعداد الصحيحة - رقمها التسلسلي في الصف X والرقم التسلسلي للصف داخل الكائن Y، وبالتالي فإن إحداثيات الزاوية اليسرى العليا للمكون هي (0، 0). . يتم تحديد عدد النقاط في كل سطر وعدد الخطوط حسب حجم المكون.

الشكل 10 - إحداثيات كائن PictureBox

في التين. يوضح الشكل 10 موقع المحاور وإحداثيات نقاط زاوية الكائن.

يتم استخدام مكون Line() لرسم محاور وأجزاء الرسم البياني المضلع للدالة.

يتلخص جوهر إنشاء الرسم البياني في حقيقة أنه يجب تقديم الوظيفة في شكل جدول (مجدول)، ثم تحديد النقاط في قالب الرسم البياني وربطها معًا.

يظهر الشكل 12 خوارزمية إنشاء الرسم البياني للوظيفة. يمكن تعديل الخوارزمية. وعلى وجه الخصوص، قد يتم الجمع بين بعض الإجراءات وقد يتم تغيير الإجراء في بعض الحالات.

دعونا نلقي نظرة على الخوارزمية بمزيد من التفصيل.

قبل تنفيذ الخوارزمية، من الضروري وصف وظيفة الروتين الفرعي لرسم الرسم البياني. وهذا ضروري لتسهيل تعديل البرنامج. إذا كنت بحاجة إلى رسم رسم بياني لوظيفة أخرى، ما عليك سوى تغيير الروتين الفرعي.

وأيضًا، قبل رسم رسم بياني، تحتاج إلى إنشاء نموذج وتحريره. يظهر في الشكل 11 مثال لتطوير نموذج. في النموذج، تحتاج إلى وضع مكونات لإدخال البيانات الأولية، ومكون لطباعة جدول، وزر أمر، وحاوية لوضع رسم بياني (PictureBox). داخل PictureBox، تحتاج إلى رسم محاور الإحداثيات باستخدام الخطوط المستقيمة ووضع التسميات لتسجيل حدود مقطع قيم وسيطة الدالة والدالة القصوى على المقطع المعني.

يتم إدخال البيانات الأولية في البرنامج المعني عن طريق الضغط على زر الأمر. في كثير من الأحيان، يتم تنفيذ إدخال البيانات باستخدام مكون TextBox.

من السهل تنفيذ إجراء جدولة الدالة في حلقة تحتوي على معلمة، نظرًا لأن عدد نقاط الرسم البياني التي سيتم حسابها معروف. قبل تنفيذ الإجراء، يجب عليك تعيين عدد صفوف الجدول.

يتم حساب عدد الأسطر بواسطة الصيغة k=n+2، حيث k هو عدد الأسطر، وn هو عدد أقسام الجدولة. يجب أن يكون عدد الأسطر أكبر من عدد المقاطع بمقدار 2، حيث من الضروري مراعاة نقطة البداية (صفر) وخط تسجيل عناوين أعمدة الصفحة.

في إجراء الجدولة نفسه، يمكنك الجمع بين إجراءين - الجدولة وحساب التطرف. يظهر هذا الحل في قائمة البرامج في الشكل 13.

تتمثل الصعوبة الرئيسية في التخطيط في الانتقال من القيمة الرياضية للدالة والوسيطة إلى إحداثيات الشاشة المستخدمة لرسم المخطط. عند حل هذه المشكلة، من الضروري مراعاة الاتجاه المعاكس للمحاور على الرسم البياني الرياضي وعلى كائن PictureBox والحاجة إلى تغيير حجم الصورة.

يتم حساب عوامل قياس المخطط باستخدام الصيغ التالية:

حيث kx هو عامل القياس على طول محور OX،

NPX – عدد وحدات البكسل الخاصة بكائن PictureBox المخصصة للتخطيط أفقيًا،

أ - القيمة الأولية لقطعة وسيطة الدالة،

ب - القيمة النهائية لقطعة وسيطة الدالة.

,

حيث Ky هو عامل القياس على طول محور OY،

NPY – عدد وحدات البكسل في كائن PictureBox المخصصة لإنشاء رسم بياني عمودي،

دقيقة - الحد الأدنى لقيمة الوظيفة،

ماكس – الحد الأقصى لقيمة الوظيفة.

تتم ترجمة الإحداثيات الرياضية للنقطة الحالية إلى إحداثيات الشاشة باستخدام الصيغ:

zx = جولة (الثور + (x(i) - أ) * kx)،

zy = Round(oy - (y(i) - Min) * ky)،

حيث zx، zy - إحداثيات الشاشة للنقطة الحالية،

ox, oy - إحداثيات نقطة تقاطع المحاور في مكون PictureBox،

x(i), y(i) – الإحداثيات الرياضية للنقطة الحالية،

kx، ky – عوامل القياس.

في صيغة حساب إحداثيات الشاشة للنقطة الحالية، يتم استخدام علامة الطرح لمراعاة الاتجاه المعاكس للمحاور (على الشاشة وعلى الرسم البياني).

تظهر قائمة البرنامج لرسم دالة في الشكل 13.

تظهر أمثلة النماذج مع نتائج البرنامج لمختلف البيانات الأولية في الشكلين 14 و 15.

الشكل 11 - مثال على تطوير النموذج

الشكل 12 - خوارزمية رسم الرسم البياني للدالة

وصف المتغيرات

Dim x() كمفرد، y() كمفرد

خاصة كمفردة

حل مسألة كوشي لنظام ODE بطريقة Runge-Kutta من الدرجة الرابعة
خوارزمية تسلسلية
صعوبة تسلسلية 4 مليون
حجم بيانات الإدخال م+3
حجم الإخراج (م+1)ن
خوارزمية متوازية
ارتفاع الشكل المتوازي المتدرج على)
عرض النموذج المتوازي المتدرج يا (م)

المؤلفون الرئيسيون للوصف: أ.أ. فورستر (1.1،1.2،1.3،1.4،1.5،1.6،1.7،1.8،1.9)، د. عليموف (1.1,1.7,1.10,2.4,2.7)

1 خصائص وهيكل الخوارزمية

1.1 وصف عام للخوارزمية

طريقة Runge-Kutta من الدرجة الرابعة هي الطريقة الأكثر شيوعًا من عائلة طريقة Runge-Kutta للخوارزميات الرقمية لحل المعادلات التفاضلية العادية وأنظمتها. تم تطوير هذه الأساليب التكرارية للحسابات التقريبية الصريحة والضمنية في حوالي عام 1900 من قبل علماء الرياضيات الألمان K. Runge وM. V. Kutta.

رسميًا، طريقة Runge-Kutta هي طريقة أويلر معدلة ومصححة، وهي عبارة عن مخططات دقة من الدرجة الثانية. هناك مخططات قياسية من الدرجة الثالثة لا يتم استخدامها على نطاق واسع.
الفكرة الرئيسية لخوارزميات Runge-Kutta هي استبدال الجانب الأيمن من المعادلة التفاضلية، التي تعتمد على الدالة المجهولة المجهولة، ببعض التقريب. إذا تمت إعادة كتابة مسألة كوشي في شكل تكامل

y(x) = y_(0) + \int\limits_(x_(0))^(x)f(t,y)dt,

ومن ثم، باستخدام الصيغ العددية المختلفة لحساب التكامل على الجانب الأيمن من المعادلة، يمكن للمرء الحصول على طرق رونج-كوتا ذات الطلبات المختلفة.

عرض عام لصيغ طرق Runge-Kutta مع خطوة الشبكة h_(n):

ص = و(ر،ص)،\،\،\، ص(t_(0)) = ذ_(0)، y_(n+1) = y_(n) + h_(n+1)\sum\limits_(i=1)^(s)b_(i)K_(n)^(i), K_(n)^(i) = f(t_(n) + c_(i)h_(n+1), y_(n) + h_(n+1)\sum\limits_(j=1)^(i) -1)أ_(ي)ك_(ن)^(ي)).

يتم تحديد طريقة Runge-Kutta المحددة من خلال مجموعة من المعاملات ب_(ط)، ج_(ي)، أ_(ي)والتي يجب أن تلبي علاقات معينة.

1.2 الوصف الرياضي للخوارزمية

1.2.1 طريقة Runge-Kutta من الدرجة الرابعة لمسألة كوشي من الدرجة الأولى DE

دعونا نفكر في مشكلة كوشي، حيث يلبي الطرف الأيمن شروط الوجود ونظريات التفرد للحل.

y" = f(x,y),\ a \leq x \leq b;\ y(a) = y^0

دعونا نضع شبكة موحدة

x_i = a + ih,\ i = 1,\dots, n,\ h = \frac(b-a)(n)

دعونا نقدم الترميز y(x_i) = y_i. نحصل على صيغة الحساب:

\begin(cases) k_1 = hf(x_i,y_i)\\ k_2 = hf(x_i + h/2,y_i + k_1/2)\\ k_3 = hf(x_i + h/2,y_i + k_2/2)\ \ k_4 = hf(x_i + h,y_i + k_3)\\ y_(i+1) = y_i + [ k_1 + 2k_2 + 2k_3 + k_4 ]/6 \\ \end(حالات)

1.2.2 طريقة Runge-Kutta من الدرجة الرابعة لمسألة كوشي لنظام DE من الدرجة الأولى

تم البحث عن حل عددي لمسألة كوشي لأنظمة ODE من الدرجة الأولى باستخدام طرق Runge-Kutta باستخدام نفس الصيغ المستخدمة في ODEs من الدرجة الأولى. على سبيل المثال، يمكن إيجاد حل بطريقة Runge-Kutta من الدرجة الرابعة إذا وضعنا:

y_i\rightarrow\bar y_i f(x_i,y_i) \rightarrow \bar f(x_i,\bar y_i) k_l \rightarrow \bar k_l \bar k_l = \begin(pmatrix) k^i_(l,1)\\ \vdots\\ k^i_(l,m)\\ \end(pmatrix)

حيث m هو البعد للنظام، ل = 1، \النقاط، 4. ونتيجة لذلك نحصل

\begin(cases) \bar k_1 = h\bar f(x_i,\bar y_i)\\ \bar k_2 = h\bar f(x_i + h/2,\bar y_i + \bar k_1/2)\\ \ bar k_3 = h\bar f(x_i + h/2,\bar y_i + \bar k_2/2)\\ \bar k_4 = h\bar f(x_i + h,\bar y_i + \bar k_3)\\ \ bar y_(i+1) = \bar y_i + [ \bar k_1 + 2\bar k_2 + 2\bar k_3 + \bar k_4 ]/6 \\ \end(حالات)

1.3 النواة الحسابية للخوارزمية

1.4 البنية الكلية للخوارزمية

تتكون البنية الكلية للخوارزمية من حساب المعاملات k_(j),\,j=1,..,4 عند العثور على قيم الوظيفة المطلوبة في كل عقدة.

1.5 مخطط تنفيذ الخوارزمية المتسلسلة

دعونا نقدم تنفيذًا محتملاً لخوارزمية متسلسلة في Matlab (يتم إعطاء الجانب الأيمن من ODE كدالة func(x,y))

1 دالة = RungeKutta4 (y0, a, b, n) 2 h = (b - a) / n; 3 x = أ : ح : ب 4 5 ص (:, 1 ) = y0 ; 6 لـ i = 2 : n 7 k1 = h * func (x (i), y (:, i)) 8 k2 = h * func (x (i) + h / 2, y (:, i) + k1 / 2 ) 9 k3 = h * func (x (i ) + h / 2 , y (:, i ) + k2 / 2 ) 10 k4 = h * func (x (i ) + h , y (:, i ) + k3 ) 11 y (:, i + 1 ) = y (:, i ) + (k1 + 2 * k2 + 2 * k3 + k4 ) / 6 12 النهاية 13 النهاية

1.6 التعقيد المتسلسل للخوارزمية

تتكون كل خطوة من حلقة الخوارزمية من 4 استدعاءات للدالة \bar f و 11 عملية ضرب و 10 عمليات إضافة. نظرًا لأن استدعاء الدالة \bar f هو العملية الأكثر تعقيدًا، فيمكن تعريف مدى تعقيد التنفيذ الخطي للخوارزمية بـ 4mn.

1.7 الرسم البياني للمعلومات

دعونا نصف الرسم البياني الخوارزمي تحليليا ورسوميا.
دعونا نعطي نظام المعادلات التفاضلية أحادية البعد. يحتوي الرسم البياني على بنية ثلاثية الأبعاد تمثل المستويات (المستويات) المترابطة. في كل مستوى i، يتم حساب التقريب i لمتجه حل النظام. هناك n من هذه المستويات في المجموع، حيث n هو عدد العقد التي يتم فيها حساب الحل التقريبي لنظام المعادلات التفاضلية. بيانات الإخراج للمستوى i هي بيانات الإدخال للمستوى i + 1.
تنقسم رؤوس الرسم البياني للمعلومات إلى نوعين:

  • النوع الأول من القمم يتوافق مع حساب إحداثيات واحدة للدالة على الجانب الأيمن من المعادلة التفاضلية عند النقاط المقدمة لمدخل الرأس. نتيجة العملية هي إحداثيات متجهة واحدة ك_(ل),\,\,\,ل = 1,2,3,4. في كل مستوى سيكون هناك 4 ملايين من هذه القمم، ويبلغ إجمالي عددها 4 ملايين.
  • النوع الثاني من القمم يتوافق مع العملية a + bc. ويمكن أيضًا تقسيم هذه القمم إلى مجموعتين. عند رؤوس المجموعة الأولى، يتم حساب التعبيرات، وسيتم توفير نتيجتها كوسائط للقمم من النوع الأول، وستكون البيانات المدخلة لهذه القمم هي نتائج إثارة رؤوس m من النوع الأول، الموجودة في نفس المستوى. في كل مستوى من رؤوس المجموعة الأولى سيكون هناك 3م، وسيكون عددهم الإجمالي 3مليون. ستكون البيانات المدخلة لرؤوس المجموعة الثانية نتيجة تشغيل قمة واحدة من النوع الأول وقمة مماثلة من النوع الثاني. ستكون نتيجة تشغيل هذا الرأس على المستوى i هي إحداثيات m لمتجه التقريب i لحل النظام. في كل مستوى سيكون هناك 4 ملايين من هذه القمم، ويبلغ إجمالي عددها 4 ملايين.

لنقدم رسمًا توضيحيًا للرسم البياني للمعلومات (الشكل 1). من أجل عدم فوضى الرسم البياني بشكل مفرط، فإننا نعتبر حالة نظام من معادلتين تفاضليتين. يشار إلى رؤوس النوع الأول باللون الأحمر، والقمم من النوع الثاني باللون الأخضر. يتم توفير البيانات المدخلة للمشكلة إلى جميع الرؤوس من النوع الأول، وكذلك إلى الرؤوس m اليسرى للمجموعة الثانية من النوع الثاني، ويتم توفير القيم الأولية للدالة المطلوبة إلى جميع الرؤوس من النوع الثاني المجموعة الأولى من النوع الثاني في المستوى الأول (انظر الفقرة). ستكون بيانات الإخراج هي نتائج إطلاق الرؤوس اليمنى للمجموعة الثانية من النوع الثاني عند كل مستوى. نتائج تشغيل القمم المتبقية هي بيانات وسيطة للخوارزمية.

الشكل 1. الرسم البياني للمعلومات لنظام من معادلتين تفاضليتين.

1.8 موارد توازي الخوارزمية

نظرًا لأنه في المخطط الحسابي الموضح أعلاه، فإن العملية الأكثر كثافة في العمالة هي حساب الجوانب اليمنى من ODE عند الحساب k_i (i = 1، \dots، 4)، فيجب إيلاء الاهتمام الرئيسي لموازاة هذه العملية. هنا سوف نستخدم نهج تحليل معادلات نظام ODE إلى أنظمة فرعية. لذلك، للتهيئة، ضع في الاعتبار مخطط تحليل البيانات التالي لعناصر المعالج المتوفرة مع الذاكرة المحلية: لكل \mu - PE (عنصر المعالج) ( \mu = 0، \dots، p-1) يتم توزيعها بواسطة المعادلات التفاضلية n/p والمتجه \bar y_0 . يتم إجراء المزيد من الحسابات وفقًا للمخطط التالي:

  1. في كل PE، يتم حساب المكونات المقابلة n/p للمتجه \bar k_1 في وقت واحد باستخدام الصيغة [ \bar k_1 ]_(\mu) = h[ \bar f(x_i, \bar y_i) ]_(\mu)
  2. لضمان مرحلة الحساب الثانية، من الضروري تجميع المتجه \bar k_1 بالكامل على كل PE. ثم يتم حساب مكونات المتجه \bar k_2 بشكل مستقل باستخدام الصيغة [ \bar k_2 ]_(\mu) = h[ \bar f(x_i + h/2,\bar y_i + \bar k_1/2)]_(\mu);
  3. يتم تجميع المتجه \bar k_2 على كل PE، ويتم حساب مكونات المتجه \bar k_3:\ [ \bar k_3 ]_(\mu) = h [\bar f(x_i + h/2,\bar y_i + \bar k_2/2)]_(\mu);
  4. يتم تجميع المتجه \bar k_3 على كل PE، ويتم حساب مكونات المتجه \bar k_4:\ [ \bar k_4 ]_(\mu) = h [\bar f(x_i + h,\bar y_i + \bar k_3)]_(\mu);
  5. يتم حساب مكونات المتجهات بالتوازي التام \bar y_(i+1):\ [\bar y_(i+1)]_(\mu) = [\bar y_(i)]_(\mu) + ([ \bar k_1 ]_(\mu ) + 2[ \bar k_2 ]_(\mu) + 2[ \bar k_3 ]_(\mu) + [ \bar k_4 ]_(\mu))/6\ويتم تجميع المتجه \bar y_(i+1) على كل PE. إذا كان من الضروري مواصلة العملية الحسابية، فسيتم افتراض i = i + 1 ويتم الانتقال إلى الخطوة 1

لاحظ أن هذه الخوارزمية لها توازي محدود، ولكن ليس توازيًا هائلًا، نظرًا لأن دورات الخوارزمية تعتمد على المعلومات.

في كل مستوى في هذه الخوارزمية، ينفذ كل PE أربع عمليات لحساب الجوانب اليمنى من ODE، وستة عشر عملية لإضافة المتجهات وضرب المتجه برقم، بالإضافة إلى نقل البيانات p-1 بين PEs الأخرى، مما يؤدي إلى إبطاء كبير جدًا أسفل تنفيذ الخوارزمية وهو عبء عند موازاة الخوارزمية. في هذه الحالة، عدد تكرارات الحلقة يساوي طول المتجه x، أي n. مما سبق يترتب على ذلك أنه عند التصنيف حسب ارتفاع PPA، يكون تعقيد الخوارزمية O(n)، وعند التصنيف حسب عرض PPA - O(m) (بشرط أن p = m).

1.9 بيانات الإدخال والإخراج للخوارزمية

بيانات الإدخال للخوارزمية هي:

  1. المتجه y^0 للبعد m;
  2. حدود الفاصل الزمني أ و ب؛
  3. معدل أخذ العينات ن؛

إجمالي حجم الإدخال: م + 3

الإخراج هو

  1. n المتجهات \bar y_i للبعد m؛
  2. المتجه \bar x للبعد n

إجمالي حجم بيانات الإخراج: (م+1)ن

1.10 خصائص الخوارزمية

دعونا ندرج الخصائص الرئيسية للخوارزمية:

2 تنفيذ البرمجيات للخوارزمية

2.1 ميزات تنفيذ الخوارزمية التسلسلية

2.2 محلية البيانات والحسابات

2.2.1 مكان تنفيذ الخوارزمية

2.2.1.1 هيكل الوصول إلى الذاكرة والتقييم النوعي للمكان
2.2.1.2 تحديد المنطقة
2.2.1.3 التحليل بناءً على اختبار Apex-Map

2.3 الطرق والميزات الممكنة للتنفيذ المتوازي للخوارزمية

2.4 قابلية التوسع للخوارزمية وتنفيذها

2.4.1 قابلية التوسع في الخوارزمية

2.4.2 قابلية التوسع في تنفيذ الخوارزمية

دعونا نجري دراسة لقابلية التوسع في التنفيذ الموازي لطريقة Runge-Kutta وفقًا للمنهجية. تم إجراء البحث على الكمبيوتر العملاق Lomonosov التابع لمجمع الكمبيوتر العملاق بجامعة موسكو. لقد درسنا التنفيذ المتوازي الخاص بنا للخوارزمية، المكتوبة بلغة C++ باستخدام معيار MPI. تم تجميع البرنامج باستخدام الإصدار 15.0 من مترجم Intel (دون تحديد مفاتيح الترجمة) باستخدام إصدار مكتبة OpenMPI 1.8.4. تم إجراء الحسابات على العقد من مجموعة T-Blade2 مع معالجات Intel Xeon 5570 Nehalem بسرعة 2.93 جيجا هرتز، مع تخصيص نواة واحدة لكل عملية MPI.

مجموعة وحدود قيم المعلمات القابلة للتغيير لبدء تنفيذ الخوارزمية:

  • عدد المعالجات في خطوات القوى اثنين حتى 64، ثم في خطوات 10؛
  • البعد للنظام.

نتيجة للتجارب، تم الحصول على النطاق التالي من كفاءة تنفيذ الخوارزمية:

  • الحد الأدنى لكفاءة التنفيذ 0.0000032%؛
  • الحد الأقصى لكفاءة التنفيذ 0.0024%.

دعونا ندرج بعض ميزات التنفيذ الموازي الذي تم اختباره:

  • للاختبار، تم استخدام الجانب الأيمن شبه العشوائي من النظام، والذي يتكون من تراكبات من الوظائف الأولية المختارة عشوائيًا؛
  • تم تعيين القيم الأولية للدالة المطلوبة بشكل عشوائي؛
  • نظرًا لأن الجانب الأيمن من النظام غير معروف في الحالة العامة، فمن غير الممكن قياس الأداء بالجيجا فلوب؛ تم استخدام عدد النقاط التي تم تقييم الوظيفة فيها كمقياس للأداء.
  • الشكل 3. كفاءة الخوارزمية المتوازية حسب حجم النظام وعدد المعالجات.

    دعونا نبني تقديرات لقابلية التوسع في التنفيذ المتوازي الذي تم اختباره لطريقة Runge-Kutta من الدرجة الرابعة:

    • حسب عدد العمليات -0.000000768. مع زيادة عدد العمليات، تنخفض الكفاءة في نطاق التغييرات المدروسة في معلمات الإطلاق، ولكن النقصان بطيء جدًا، وهو ما يرتبط بالحد الأقصى المنخفض للغاية لكفاءة التنفيذ المتوازي للخوارزمية. يتم تفسير انخفاض الكفاءة مع زيادة عدد العمليات للأنظمة ذات البعد الصغير بحقيقة أنه إذا تجاوز عدد العمليات أبعاد النظام، فلن يشارك بعضها في الحسابات. بالنسبة للأنظمة ذات البعد الكبير، تؤثر الزيادة في العمليات سلبا على الكفاءة، حيث يزيد عدد عمليات النقل بينهما، مما يبطئ العمل بشكل كبير.
    • حسب حجم المهمة 0.000000763. مع زيادة أبعاد النظام مع عدد ثابت من العمليات، تزداد كفاءة التغييرات في معلمات الإطلاق في المنطقة قيد النظر، ويتم ملاحظة الزيادة في المنطقة قيد النظر بأكملها تقريبًا. يتم تفسير ذلك بحقيقة أنه مع عدد ثابت من العمليات، تؤدي الزيادة في أبعاد النظام إلى زيادة الحمل الحسابي على كل معالج، ونتيجة لذلك يزداد الوقت الحسابي بالنسبة إلى وقت الخمول. مع زيادة عدد العمليات، ينخفض ​​معدل الزيادة في الكفاءة مع زيادة حجم المهمة، حيث يصبح عدد عمليات النقل كبيرًا جدًا.
    • في اتجاهين 0.00000000174. ومع الزيادة المتزامنة في عدد العمليات وحجم النظام، تزداد الكفاءة. معدل الزيادة في الكفاءة منخفض للغاية، وهو ما يفسره ارتفاع التكاليف العامة.

    2.5 الخصائص الديناميكية وكفاءة تنفيذ الخوارزمية

    2.6 استنتاجات لفصول الهندسة المعمارية

    2.7 التطبيقات الحالية للخوارزمية

    طريقة الترتيب الرابع هي الأكثر استخدامًا بين جميع مخططات Runge-Kutta، لذلك هناك العديد من التطبيقات البرمجية المتسلسلة لها، سواء التجارية أو المجانية. أحد أشهر تطبيقات البرامج هو ode45 في بيئة MATLAB، والذي يحتوي على عدد كبير من الخيارات لتخصيص الحسابات الرقمية، مما يجعله مناسبًا تمامًا للاستخدام. كما يتم أيضًا تطبيق طريقة Runge-Kutta من الدرجة الرابعة في مكتبة PETSc الموازية، والتي يتم توزيعها مجانًا. على الرغم من أن معظم الأساليب في هذه المكتبة يتم تنفيذها باستخدام خوارزميات متوازية، إلا أن طريقة Runge-Kutta يتم تنفيذها بشكل تسلسلي. من الصعب جدًا العثور على تطبيق موازٍ للطريقة عند النظر في مشكلة كوشي العامة، ولكن هناك العديد من الأعمال العلمية التي تصف التطبيقات المتوازية لطريقة Runge-Kutta لفئات معينة من الأنظمة، على سبيل المثال، في أعمال A.V. يصف ستارتشينكو التنفيذ المتوازي للأنظمة الخطية.

طريقة Runge-Kutta الكلاسيكية 4 أوامر

الموضوع 6.3. طريقة رونج-كوتجي

طرق رونج-كوتا- عائلة مهمة من الخوارزميات العددية لحل (أنظمة) المعادلات التفاضلية العادية. تم تطوير هذه الأساليب التكرارية للحسابات التقريبية الصريحة والضمنية في حوالي عام 1900 من قبل علماء الرياضيات الألمان K. Runge وM. V. Kutta.

رسميًا، طرق رونج-كوتا هي طريقة أويلر معدلة ومصححة، وهي تمثل مخططات دقة من الدرجة الثانية. هناك مخططات قياسية من الدرجة الثالثة لا يتم استخدامها على نطاق واسع. الأكثر استخدامًا وتنفيذًا في الحزم الرياضية المختلفة (Maple، MathCAD، Maxima) هو مخطط الترتيب الرابع القياسي (انظر أدناه). في بعض الأحيان، عند إجراء العمليات الحسابية بدقة متزايدة، يتم استخدام مخططات الترتيب الخامس والسادس.

إن طريقة Runge-Kutta من الدرجة الرابعة منتشرة على نطاق واسع لدرجة أنها غالبًا ما تسمى ببساطة طريقة Runge-Kutta.

خذ بعين الاعتبار مشكلة كوشي. ثم يتم حساب القيمة عند النقطة التالية باستخدام الصيغة:

حجم خطوة الشبكة هو .

هذه الطريقة لها الترتيب الرابع، أي. الخطأ في كل خطوة هو O(h5)، والخطأ الإجمالي في فترة التكامل النهائية هو O(h4).

عائلة طرق Runge-Kutta المباشرة هي تعميم لطريقة Runge-Kutta من الدرجة الرابعة. يتم تقديمه بواسطة الصيغ

يتم تحديد الطريقة المحددة من خلال الرقم s والمعاملات bi وaij وci. غالبًا ما يتم تنظيم هذه المعاملات في جدول

دعونا نعطي معادلة تفاضلية من الدرجة الأولى مع الشرط الأولي y(x 0)=y 0. دعونا نختار الخطوة h ونقدم الترميز:

x i = x 0 + ih و y i = y(x i)، حيث i = 0, 1, 2, ... .

على غرار الطريقة الموصوفة أعلاه، يتم إجراء الحل

المعادلة التفاضلية. الفرق هو أن الخطوة مقسمة إلى 4 أجزاء.

وفقًا لطريقة Runge-Kutta من الترتيب الرابع، يتم تحديد القيم المتعاقبة y i للدالة المطلوبة y بواسطة الصيغة:

ذ أنا+1 = ذ أنا +؟ص أنا

حيث أنا = 0، 1، 2 ...

ص=(ك1+2*ك2+2*ك3+ك4)/6

يتم حساب الأرقام k1، k2، k3، k4 في كل خطوة باستخدام الصيغ:

k1 = ح*و(س ط، ص ط)

k2 = و (x i +h/2, y i +k1 /2)*h

k3 = F(x i +h/2, y i +k2 /2)*h

k4 = F(x i +h, y i +k3)*h

هذه طريقة واضحة من أربع خطوات مع 4 أوامر من الدقة.

يظهر في الشكل 6 مخطط تدفق الإجراء الخاص بحل المعادلة التفاضلية باستخدام طريقة Runge-Kutta.

F(x, y) - دالة معينة - يجب وصفها بشكل منفصل.

معلمات الإدخال:

X0، XK - الأولي والنهائي

القيم المستقلة

عامل؛

Y0 - القيمة y 0 من الحالة الأولية

ن - عدد أجزاء القسم؛

معلمات الإخراج:

Y - مجموعة قيم الحل المطلوب

في العقد الشبكة.

خوارزمية لحل المشكلة

الحل في MathCAD

قائمة البرامج بلغة Visual Basic

Dim xr(), yr(), xe(), ye(), xo(), yo() كمفردة

خاص x0، y0، h، xk، k1، k2، c، k3، k4، yd كفردي

خاص n، i كعدد صحيح

الوظيفة العامة f (ByVal a، ByVal b) كمفردة

و = -(أ - ب) / أ

الزر الفرعي الخاص 1_Click (مرسل ByVal كـ System.Object، ByVal e كـ System.EventArgs) يعالج Button1.Click

x0 = فال (TextBox1.Text)

xk = فال (TextBox2.Text)

y0 = فال (TextBox4.Text)

ح = فال (TextBox3.Text)

ن = (xk - x0) / ح

ج = y0 / x0 + Math.Log(x0)

DataGridView1.ColumnCount = 4

DataGridView1.RowCount = n + 2

DataGridView1.Item(0, 0).Value = "x"!}

DataGridView1.Item(1, 0).Value = "عام"!}

DataGridView1.Item(2, 0).Value = "Euler M"!}

DataGridView1.Item(3, 0).Value = "Runge Kutt"!}

لأني = 0 إلى ن - 1

xe(i) = Math.Round((xe(0) + i * h), 2)

انتم(i + 1) = انتم(i) + h * f(xe(i) + h / 2, ye(i) + h / 2 * f(xe(i), انتم(i)))

DataGridView1.Item(2, 1).القيمة = نعم(0)

DataGridView1.Item(0, 1).القيمة = xe(0)

DataGridView1.Item(0, i + 1).Value = xe(i)

DataGridView1.Item(2, i + 1).Value = Str(ye(i))

لأني = 0 إلى ن - 1

xr(i) = Math.Round((xe(0) + i * h), 2)

k1 = ح * و(xr(i)، سنة(i))

k2 = ح * و(xr(i) + ح / 2، سنة(i) + k1 / 2)

k3 = ح * و(xr(i) + ح / 2، سنة(i) + k2 / 2)

k4 = ح * و(xr(i) + ح، سنة(i) + k3)

ياردة = (k1 + 2 * k2 + 2 * k3 + k4) / 6

سنة (ط + 1) = سنة (أنا) + ياردة

DataGridView1.Item(3, 1).القيمة = سنة(0)

DataGridView1.Item(3, i + 1).Value = Str(yr(i))

لأني = 0 إلى ن - 1

xo(i) = Math.Round((xe(0) + i * h), 2)

yo(i) = xo(i) * (ج - Math.Log(xo(i)))

DataGridView1.Item(1, 1).القيمة = يو(0)

DataGridView1.Item(1, i + 1).Value = Str(yo(i))

Chart1.Series.Add("الحل العام")

Chart1.Series("الحل العام").ChartType = SeriesChartType.Line

لأني = 0 إلى ن - 1

Chart1.Series("الحل العام").Points.AddXY(xo(i), yo(i))

Chart1.Series("الحل العام").ChartArea = "ChartArea1"

Chart1.Series.Add("أويلر م")

Chart1.Series("أويلر M").ChartType = SeriesChartType.Point

لأني = 0 إلى ن - 1

Chart1.Series("أويلر M").Points.AddXY(xe(i), ye(i))

Chart1.Series("أويلر م").ChartArea = "ChartArea1"

Chart1.Series("Euler M").اللون = اللون.أزرق

Chart1.Series.Add("رونجي كوت")

Chart1.Series("Runge Kutt").ChartType = SeriesChartType.Line

لأني = 0 إلى ن - 1

Chart1.Series("رونجي كوت").Points.AddXY(xr(i), yr(i))

Chart1.Series("رونجي كوت").ChartArea = "ChartArea1"

Chart1.Series("Runge Kutt").اللون = اللون.أخضر

وعلوم الكمبيوتر

معهد الأورال التقني للاتصالات والمعلوماتية

قسم الفيزياء والرياضيات التطبيقية وعلوم الحاسوب.

عمل الدورة

في علوم الكمبيوتر:

تصور الطرق العددية.

حل المعادلات التفاضلية العادية.

أكمله: أرابوف

غرام. مي-71

تم الفحص بواسطة: Minina E. E.

ايكاترينبرج

2008

مقدمة ……………………………………………………………………….3

1. بيان المشكلة ………………………………………….4

2. وصف طرق الحل …………………………………………..5

2. 1. جوهر المهمة ……………………………………………….5

2. 2. المعنى الهندسي للمشكلة………………………………….5

2. 3. الطرق العددية لحل مسألة كوشي ...........................6

2. 4. طريقة أويلر ………………………………………………….6

2. 5. طريقة رونج-كوتا من الدرجة الرابعة.................................8

2. 6. حل المشكلة باستخدام طرق أويلر و

رونج-كوتا الترتيب الرابع ………………………………….9

2. 6. 1. طريقة أويلر ………………………………….9

2. 6. 2. طريقة رونج-كوتا من الدرجة الرابعة .......................... 10

3. خوارزمية حل المشكلة ...........................11

3. 1. خوارزميات الإجراءات الفرعية ........................................... 11

3. 1. 1. روتين فرعي لطريقة أويلر...........................11

3. 1. 2. الروتين الفرعي لطريقة Runge-Kutta من الدرجة الرابعة ............12

3. 1. 3. روتين الحل العام ............ …………………1 3

3. 2. خوارزمية الوظيفة ………………………………………… 1 3

3. 3. خوارزمية البرنامج ……………………………………………………………………………………………………… 14

4. نموذج البرنامج.................................................................................17

5. قائمة البرامج................................................................................................. 18

6. حل المشكلة في MathCad …………………………………………..20

الخلاصة ………………………………………………… 22


مقدمة

العديد من القوانين الفيزيائية التي تحكم بعض الظواهر مكتوبة على شكل معادلة رياضية تعبر عن علاقة معينة بين كميات معينة.ترجع الأهمية الكبيرة التي تتمتع بها المعادلات التفاضلية للرياضيات وخاصة لتطبيقاتها إلى حقيقة أن دراسة العديد من المشكلات الفيزيائية والتقنية تتلخص في حل مثل هذه المعادلات. تلعب المعادلات التفاضلية أيضًا دورًا مهمًا في العلوم الأخرى، مثل علم الأحياء والاقتصاد والهندسة الكهربائية؛ في الواقع، فإنها تنشأ حيثما تكون هناك حاجة إلى وصف كمي (عددي) للظواهر. ولذلك، فإن حل المعادلات التفاضلية سيكون دائمًا مهمة ضرورية وعاجلة.

الغرض من هذا المقرر هو حل المعادلات التفاضلية باستخدام طريقتين عدديتين: طريقة أويلر وطريقة رونج-كوتا من الدرجة الرابعة من الدقة.

لتحقيق الهدف، حددت لنفسي المهام التالية:

  1. أكتب برنامج لحل هذه المعادلة التفاضلية باستخدام طريقتين عدديتين في البرنامجالبصرية الأساسية.
  2. التحقق من الحل باستخدام التطبيق MathCad.
  3. قارن النتائج التي تم الحصول عليها بطرق مختلفة مع الحل العام.


1. بيان المشكلة

حل طرق أويلر وطرق أويلر المعدلة لحل مشكلة كوشي لمعادلة تفاضلية من الدرجة الأولى على الفترة [ X0؛ X ك ] مع الخطوة ح والشرط الأولي:ص (س 0 ) = ص 0 .

وينبغي الحصول على الإجابة في شكل جدول النتائج:

نعم (1)

نعم (2)

ص 0 (1)

ص 0 (2)

ص(× 0)

ص 1 (1)

ص 1 (2)

ص(× 1)

واي ك (1)

واي ك (2)

ص (Xك)

حيث ص (1) ، ص (2) الحلول التي تم الحصول عليها بطرق عددية مختلفة،واي تي الحل الدقيق للمعادلة التفاضلية.

من الممكن عرض نتائج الحل ليس على شكل جدول، بل على شكل قوائم.

تصور بيانات الجدول في نموذج في شكل رسوم بيانية.

قبل حساب العمود الأخير من جدول النتائج، من الضروري طرح قيمة المعامل من الشروط الأوليةج ، تستخدم في الحل العام.

المعادلة التفاضلية

قرار مشترك

ص = - س·ص/(س+1)

1 ,2

0, 1

ص=ج (س+1) إكسب(-س)


2. وصف طرق الحل

2. 1. جوهر المشكلة

لحل معادلة تفاضلية عادية لا بد من معرفة قيم المتغير التابع و (أو) مشتقاته لقيم معينة للمتغير المستقل. إذا تم تحديد هذه الشروط الإضافية لقيمة واحدة من المتغير المستقل، فإن مثل هذه المشكلة تسمى مشكلة الشروط الأولية، أو مشكلة كوشي. في كثير من الأحيان في مشكلة كوشي، يلعب الزمن دور متغير مستقل.

ويمكن صياغة مسألة كوشي على النحو التالي:

دع المعادلة التفاضلية والحالة الأولية تعطىص (س 0 ) = ص 0 . نحن بحاجة إلى العثور على الدالة y(س ) ، محققًا كلاً من المعادلة المحددة والشرط الأولي.

يتم تقليل الحل العددي لمشكلة كوشي إلى جدولة الوظيفة المطلوبة.

يسمى الرسم البياني لحل المعادلة التفاضلية بالمنحنى التكاملي.

2. 2. المعنى الهندسي للمشكلة

ص = و (س، ص ) - ظل زاوية ميل المماس للرسم البياني للحل عند النقطة (x، y) إلى المحور 0X، - المعامل الزاوي (الشكل 1).

الشكل 1. المعنى الهندسي لمشكلة كوشي.

وجود الحل:

إذا كان الجانب الأيمنو(س، ص ) مستمر في بعض المناطقر ، والتي تحددها عدم المساواة

| س - س 0 |< а; | y - y 0 | < b ,

ثم هناك حل واحد على الأقل y = y(x)، محدد في الحي |x x 0 | < h , где h - رقم موجب، عدد إيجابي.

هذا الحل فريد من نوعه إذار تم استيفاء شرط Lipschitz

| و (س، ص)- و (س، ص)| ≥ ن | ص - ص |(س، ص)،

حيث N هو ثابت ما (ثابت Lipschitz)، ويعتمد، في الحالة العامة، على a وب. إذا و(س ، ذ) لديه مشتق محدود

و ذ (س، ذ) في ر ثم يمكننا وضع N = max |و ذ (س، ذ)| مع (x، y) التي تنتمي إلى R.

2. 3. الطرق العددية لحل مسألة كوشي

عند استخدام الطرق العددية، المقطع [x 0، اكس ] - مناطق التغيير المستمر للوسيطة x بواسطة مجموعة. تتكون من عدد محدود من النقاط x 0 < х 1 < ... < x n = Х - сеткой.

في هذه الحالة × ط تسمى العقد الشبكة.

تستخدم العديد من الطرق شبكات موحدة ذات مسافات متباعدة:

مشكلة كوشي، تم تعريفها مسبقًا على الفترة المستمرة [x 0، اكس ]، يتم استبداله بنظيره المنفصل - نظام من المعادلات، يمكن حله من خلال إيجاد القيم على التواليذ 1 , ص 2 ,…, ذ ن - القيم التقريبية للوظيفة في عقد الشبكة.

يستخدم الحل العددي لمسألة كوشي على نطاق واسع في مختلف مجالات العلوم والتكنولوجيا، وعدد الطرق التي تم تطويرها لها كبير جدًا. يمكن تقسيم هذه الأساليب إلى المجموعات التالية.

طرق من خطوة واحدة للعثور على النقطة التالية
منحنى ص =و(س ) مطلوب معلومات حول خطوة سابقة واحدة فقط.
طريقة أويلر وطرق رونج-كوتا هي خطوة واحدة.

طرق التنبؤ والتصحيح (متعددة الخطوات)، حيث يتم العثور على النقطة التالية من المنحنى y =و(س ) مطلوب معلومات عن أكثر من نقطة من النقاط السابقة. للحصول على قيمة عددية دقيقة بما فيه الكفاية، غالبا ما يتم استخدام التكرار. وتشمل هذه الأساليب أساليب ميلن وآدامز-باشفورث وهامينغ.

الطرق الصريحة التي لا تعتمد عليها الدالة Фص ن +1 .

الطرق الضمنية التي تعتمد عليها الدالة Фص ن +1 .

2.4 طريقة أويلر.

تسمى هذه الطريقة أحيانًا طريقة Runge-Kutta للدقة من الدرجة الأولى.

هذه الطريقة هي خطوة واحدة. يتم جدولة الوظيفة في كل نقطة على حدة. لحساب قيمة دالة في العقدة التالية، يجب عليك استخدام قيمة الدالة في عقدة واحدة سابقة.

دعونا نعطي معادلة تفاضلية من الدرجة الأولى:

ص = و(س، ص)

مع الحالة الأولية

ص (س 0 ) = ص 0

دعونا نختار الخطوة ح وإدخال التدوين:

x i = x 0 + ih و y i = y (x i )، حيث i = 0، 1، 2، ...،

س ط - العقد الشبكة،

ذ ط - قيمة الدالة المتكاملة في العقد.

وتظهر الرسوم التوضيحية للحل في الشكل 2.

لنرسم خطًا مستقيمًا AB عبر النقطة (س ط، ذ ط ) في الزاوية α. حيثتيراغرام α = و (س ط، ص ط )

وفقًا للمعنى الهندسي للمسألة، يكون الخط المستقيم AB مماسًا للدالة التكاملية. دعونا نستبدل نقطة الدالة التكاملية بنقطة تقع على المماس AB.

ثم y i +1 = y i + Δy

دعونا نساوي الجانبين الأيمنتيراغرام α = و (س ط، ص ط) و. نحن نحصل

ومن ثم Δ у = h ∙ f (x i، y i).

دعونا نستبدل الصيغة في هذا التعبيرص ط +1 = ص ط + Δy ، ثم قم بتحويله. ونتيجة لذلك، نحصل على صيغة لحساب النقطة التالية من وظيفة التكامل:

الشكل 2. طريقة أويلر.

يتضح من الصيغة أنه لحساب كل نقطة لاحقة من الدالة التكاملية، من الضروري معرفة قيمة نقطة واحدة سابقة فقط. وبالتالي، بمعرفة الشروط الأولية، من الممكن بناء منحنى متكامل على فترة زمنية معينة.

يظهر الشكل 3 مخططًا تدفقيًا لإجراءات حل المعادلة التفاضلية باستخدام طريقة أويلر.

و(خ ، ذ) - يجب أن تكون الوظيفة المحددة

سيتم وصفها بشكل منفصل.

معلمات الإدخال:

اكس0، اكس كيه البداية والنهاية

قيم متغيرة مستقلة؛

ص 0 القيمة ص 0 من الحالة الأولية

ص (س 0 ) = ص 0 ;

معلمات الإخراج:

Y - مجموعة قيم الحل المطلوب

في العقد الشبكة.

الشكل 3. مخطط تدفق الإجراء لحل المعادلة التفاضلية باستخدام طريقة أويلر.

تعد طريقة أويلر واحدة من أبسط الطرق لحل المعادلات التفاضلية العادية عدديًا. لكن عيبها الكبير هو الخطأ الكبير في الحساب. في الشكل 2، خطأ الحساب لانا اذهب يتم الإشارة إلى الخطوة بواسطة ε. ومع كل خطوة يزداد الخطأ الحسابي.

2.5 طريقة Runge-Kutta 4 أوامر

دع معادلة تفاضلية من الدرجة الأولى تعطى مع الشرط الأولي y (x 0 )= y 0. اختر الخطوة h وإدخال التدوين:

x i = x 0 + ih و y i = y (x i )، حيث i = 0, 1, 2, ... .

على غرار الطريقة الموصوفة أعلاه، يتم الحل

المعادلة التفاضلية. الفرق هو أن الخطوة مقسمة إلى 4 أجزاء.

وفقا لطريقة Runge-Kutta من الدرجة الرابعة، القيم المتعاقبة y i للوظيفة المطلوبة y يتم تحديدها بواسطة الصيغة:

y i+1 = y i +∆y i حيث i = 0, 1, 2 ...

∆y=(k1+2*k2+2*k3+k4)/6

أ الأرقام ك 1، ك 2، ك 3، ك يتم حساب 4 في كل خطوة باستخدام الصيغ:

ك 1 = ح * و (س ط , ص ط )

k2 = و (x i +h/2, y i +k1 /2)*h

k3 = F(x i +h/2, y i +k2 /2)*h

k4 = F(x i +h, y i +k3) * h

هذه طريقة واضحة من أربع خطوات مع 4 أوامر من الدقة.

يظهر في الشكل 6 مخطط تدفق الإجراء الخاص بحل المعادلة التفاضلية باستخدام طريقة Runge-Kutta.

و(خ ، ذ) - وظيفة معينة - يجب

سيتم وصفها بشكل منفصل.

معلمات الإدخال:
×0، × ك - الابتدائي والنهائي

القيم المستقلة

عامل؛

ص 0 القيمة ص 0 من الحالة الأولية

ص (س 0 )= ذ 0 ;

ن - عدد أجزاء القسم؛

معلمات الإخراج:

ي - مجموعة قيم الحل المطلوب

في العقد الشبكة.

2. 6. حل المشكلة التي تطرحها طرق أويلر ورونجي-كوتا من الدرجة الرابعة

2. 6. 1. طريقة أويلر

1. بناء محاور الإحداثيات.

2. مارك أ (1،2؛ 1) النقطة الأولى من المنحنى التكاملي؛

4. نحن نبني الظلل 0 عند النقطة A عند الزاوية α 0 ;

5. أوجد x 1 باستخدام الصيغة: x i = x 0 + ih، حيث h خطوة التكامل

× 1 = 1.2 + 1 0.1 = 1.3

6. نقوم بتنفيذ مباشرس = س 1 = 0.1 إلى التقاطع مع الخطل 0، ضع علامة على النقطة B (x 1؛ y 1)؛

7. ابحث عن النقطة y B:

من المثلث الأيمن ABC,

Δ ص = ص 1 ذ 0 ,

Δس = س 1 × 0 = ح،

f (x 0 ; y 0 ) = (y 1 y 0 )/ h =>

y 1 = y 0 + h (f (x 0 ; y 0 )) = 1 + 0.1 f (1.2;1) = 1-0.545454 = 0.945

ولذلك النقطةب له إحداثيات (1.3؛ 0.945).

2. 6 .2. طريقة Runge-Kutta 4 أوامر

1. بناء محاور الإحداثيات.

2. ضع علامة A(1,2; 1) على النقطة الأولى للمنحنى التكاملي؛

3. نحن نبحث عن زاوية ميل المماس للرسم البياني عند هذه النقطةأ:

4. نحن نبني الظلل 0 عند النقطة A عند الزاوية α 0 ;

5. أوجد x 1 باستخدام الصيغة: x i = x 0 + ih

س 1 = 1,2 + 1 0, 1 = 1, 3;

  1. نجدها باستخدام الصيغ:

ك1=0.1·و(1.2;1)=0.1·(-0.545454)= - 0.0545

k2=0.1· f(1.2+0.1/2;1-0.0545/2)= -0.0540

K3=0.1· f(1.2+0.1/2;1-0.0540/2)= - 0.0540

K4=0.1· f(1.2+0.1;1-0.0540)= - 0.0535

∆ص 1 =(- 0.0545+2·(-0.0540)+2·(-0.0540) - 0.0535)/6= - 0.054

∆ ص 2 =1- 0.054=0.946

ولذلك، فإن النقطة التالية على الرسم البياني للحل لها إحداثيات (1.3؛ 0.946).

3. خوارزمية حل المشكلة

3. 1. خوارزميات الروتين الفرعي

3.1.1 روتين فرعي لطريقة أويلر

3.1.2 الروتين الفرعي لطريقة Runge-Kutta 4 أوامر

3. 1. 3. روتين الحل العام

3. 2. خوارزمية الوظيفة


3. 3. خوارزمية البرنامج

4. استمارة البرنامج

  1. قائمة البرنامج

Dim j() كمفردة

خافت x () كمفردة

خافت y() كمفردة

خافت o () كمفردة

خاص n، i كعدد صحيح

خاص xk، x0، kx، ky كفرد

خاص ك، ك1، ك2، ك3، ك4 كفرد

خاص h، max، min، y0 كفردي

وظيفة خاصة f(a, b كمفردة) كمفردة

و = -أ * ب / (أ + 1)

وظيفة النهاية

وظيفة خاصة f1(x كمفردة) كمفردة

وظيفة النهاية

ايلر الفرعية الخاصة()

إعادة ديم س (ن)

رديم ي (ن)

ي(0) = ذ0

لأني = 0 إلى ن

س(ط) = س0 + ح * ط

بعدها انا

لأني = 0 إلى ن - 1

بعدها انا

نهاية الفرعية

رونج فرعي خاص ()

إعادة ديم ص (ن)

ص(0) = ذ0

لأني = 0 إلى ن

س(ط) = س0 + ح * ط

بعدها انا

لأني = 0 إلى ن - 1

k1 = ح * و(س(ط)، ص(ط))

ص(ط + 1) = ص(ط) + ك

MSFlexGrid1.TextMatrix(1, 2) = Str(y0)

بعدها انا

نهاية الفرعية

Obhee الفرعية الخاصة ()

إعادة ديم س (ن)

لأني = 0 إلى ن

س(0) = ص0

س(ط) = س0 + ح * ط

س(ط) = f1(س(ط))

بعدها انا

نهاية الفرعية

أمر فرعي خاص1_Click()

x0 = فال(نص1.نص)

xk = فال(Text2.Text)

ح = فال (نص 4. نص)

y0 = فال(Text3.Text)

ن = (xk - x0) / ح

Label6.Caption = Str(x0)

Label5.Caption = Str(xk)

MSFlexGrid1.Rows = n + 2

MSFlexGrid1.TextMatrix(0, 1) = "صحيح"

MSFlexGrid1.TextMatrix(0, 2) = "الاستجابة"

MSFlexGrid1.TextMatrix(0, 3) = "نص آخر"

ايلر

رونج

أوبهي

الحد الأقصى = ص0

دقيقة = ص0

لأني = 0 إلى ن

إذا j (i)> الحد الأقصى ثم

الحد الأقصى = ي(ط)

إنهاء إذا

إذا ي(ط)< min Then

دقيقة = ي (ط)

إنهاء إذا

إذا y(i)> الحد الأقصى ثم

الحد الأقصى = ص (ط)

إنهاء إذا

إذا ذ (أنا)< min Then

دقيقة = ص (ط)

إنهاء إذا

إذا o(i)> الحد الأقصى ثم

الحد الأقصى = س (ط)

إنهاء إذا

إذا س (ط)< min Then

دقيقة = س (ط)

إنهاء إذا

بعدها انا

Label4.Caption = Str(max)

Label7.Caption = Str(min)

Picture1.Cls

لأني = 1 إلى ن - 1

X1 = 720 + جولة(kx * (x(i - 1) - x0))

X2 = 720 + جولة(kx * (x(i) - x0))

X1 = 720 + جولة(kx * (x(i - 1) - x0))

X2 = 720 + جولة(kx * (x(i) - x0))

بعدها انا

نهاية الفرعية

أمر فرعي خاص2_Click()

نهاية الفرعية

  1. حل المشكلة في MathCad

أويلر

رونج كوت

عام

خاتمة

من بين الطريقتين (أويلر ورونجي-كوتا)، وفقًا للنتائج التي تم الحصول عليها، تبين أن طريقة رونج-كوتا أكثر دقة (مقارنة بالحل العام). يتم تفسير ذلك بحقيقة أنه، على عكس طريقة أويلر، في طريقة Runge-Kutta، لا يتم تقسيم الخطوة إلى 4 أجزاء، ونتيجة لذلك يصبح خطأ الطريقة أصغر.

عند الانتهاء من أعمال الدورة، أكملت جميع المهام المكلفة بها: قمت بكتابة برنامج لحل هذه المعادلة التفاضلية باستخدام طريقتين عدديتين في البرنامجالبصرية الأساسية ، تحقق من الحل باستخدام التطبيق MathCad ، مقارنة النتائج التي تم الحصول عليها بطرق مختلفة مع الحل العام. أعتقد أنني قد حققت الهدف من هذه الدورة التدريبية بالكامل.


تيراغرام (α) = و (س، ص)

Y i+1 = Y i + h ∙ F(x, Y i )

س = X0 + أنا ∙ ح

ط = 0، …، ن - 1

ح = (Xk X0)/N

نهاية

MSFlexGrid1.TextMatrix(i + 2, 2) = Str(y(i + 1))

MSFlexGrid1.TextMatrix(1, 2) = Str(ص (0))

ك س = (6600 - 720) / (س ك - س 0)

كي = (6600 - 1120) / (الحد الأقصى - الحد الأدنى)

Label4.Caption = Str(max)

Label7.Caption = Str(min)

س(ط)=دقيقة

س (ط)

MSFlexGrid1.TextMatrix(i + 1, 3) = Str(o(i))

ص (ط) = الحد الأقصى

آيلر (X0، Xk، Y0، N، Y)

س ط +1

× ط

ص ط+1

ص = ص (س)

ك2=ح*و(س+ح/2، ص ط +ك1/2)

K1=ح*F(س،ص ط)

ص (ط) > الحد الأقصى

ي(ط)=دقيقة

س = X0 + أنا ∙ ح

ط = 0، …، ن-1

ح = (Xk X0)/N

Rynge4(X0، Xk، Y0، N، Y)

الأمر2

س(ط)=الحد الأقصى

س (ط) > الحد الأقصى

ك=(ك1+2*ك2+2*ك3+ك4)/6

k4= ح*F(x+h، Y i +k3)

k3= ح*F(س+ح/2، ص ط +ك2/2)

ص ط+1 = ص ط +ك

k1 = ح * و(س(ط)، ص(ط))

k2 = ح * و(x(i) + ح / 2، ص(i) + k1 / 2)

k3 = ح * و(x(i) + ح / 2، ص(i) + k2 / 2)

k4 = ح * و(س(ط) + ح، ص(ط) + ك3)

ك = (ك1 + 2 * ك2 + 2 * ك3 + ك4) / 6

ص(ط + 1) = ص(ط) + ك

ط = 0، …، ن-1

س(ط) = س0 + ح * ط

ط = 0، …، ن

إعادة ديم ص (ن)

ز(0) = ص0

يبدأ

ص0،X0،XK،ح

ن = (xk - x0) / ح

الحد الأقصى = ص0

دقيقة = ص0

ي(ط) > الحد الأقصى

أنا = 0، ... ن

ايلر

رونج

اوبشي

MSFlexGrid1.Rows = n + 2

MSFlexGrid1.TextMatrix(0, 0) = "x"

MSFlexGrid1.TextMatrix(0, 1) = "أويلر "

MSFlexGrid1.TextMatrix(0, 2) = "رونج- كوتا"

MSFlexGrid1.TextMatrix(0, 3) = "قرار مشترك"

Label6.Caption = Str(x0)

Label5.Caption = Str(xk)

ي(ط)=الحد الأقصى

X1 = 720 + جولة(kx * (x(i - 1) - x0))

X2 = 720 + جولة(kx * (x(i) - x0))

Y1 = 6600 - جولة(ky * (o(i - 1) - min))

Y2 = 6600 - جولة(ky * (o(i) - min))

نهاية

ي (ط)

X1 = 720 + جولة(kx * (x(i - 1) - x0))

X2 = 720 + جولة(kx * (x(i) - x0))

Y1 = 6600 - جولة(ky * (y(i - 1) - min))

Y2 = 6600 - جولة(ky * (y(i) - min))

X1 = 720 + جولة(kx * (x(i - 1) - x0))

X2 = 720 + جولة(kx * (x(i) - x0))

Y1 = 6600 - جولة(ky * (j(i - 1) - min))

Y2 = 6600 - جولة(ky * (j(i) - min))

أنا =1 ، …، ن-1

نهاية

f1 = y0 / ((x0 + 1) * إكسب(-x0)) * (س + 1) * إكسب(-x)

f1(خ)

MSFlexGrid1.TextMatrix(1, 0) = Str(x0)

MSFlexGrid1.TextMatrix(i + 2, 0) = Str(x(i + 1))

MSFlexGrid1.TextMatrix(i + 2, 1) = Str(j(i + 1))

MSFlexGrid1.TextMatrix(1, 1) = Str(j(0))

ي(i + 1) = ي(i) + ح * و(x(i)، ي(i))

ط = 0، …، ن-1

س(ط) = س0 + ح * ط

ط = 0، …، ن

إعادة ديم س (ن)

رديم ي (ن)

ي(0) = ذ0

ايلر

نهاية

س(ط) = س0 + ح * ط

س(ط) =f1(س (ط))

ط = 0، …، ن

إعادة ديمس(ن)

س(0) = ص0

أوبهي

نهاية

MSFlexGrid16

الصورة 1

و = -أ * ب / (أ + 1)

و(أ،ب)

رونج

لابي71

النص2

النص 1

لابي41

لابي31

التسمية1

النص3

لابي21

التسمية6

النص4

الأمر1

التسمية4

لابي51

لابي91

التسمية10

نهاية

Picture1.Line (X1، Y1)-(X2، Y2)، RGB (0، 200، 0)

Picture1.Line (X1، Y1)-(X2، Y2)، RGB (500، 70، 90)

Picture1.Line (X1، Y1)-(X2، Y2)، RGB (400، 100، 12)

ص (ط)

ذ(ط)=دقيقة