معادلات دیفرانسیل و انتگرال در متلب

ماشین حساب

خرداد ۱۶, ۱۳۹۵Kamaliافزونه نرم افزار متلب

مقدمه

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

برای مثال، وقتی حاصل دستور (sin(π/۲ را در فضای عددی متلب بکار ببریم، جواب حاصل عددی نزدیک به صفر (مخالف صفر) خواهد بود. درصورتیکه در واقعیت این مقدار، دقیقاً صفر است. حال برای اینکه مقدار واقعی را بدست آوریم، باید از ویژگی محاسبات نمادین موجود در متلب استفاده کنیم.

در حقیقت در محاسبات عددی، به خاطر خطای ذخیره‌سازی و محدودیت حافظه، هر محاسبه‌ای که انجام می‌شود با اندکی تقریب در اختیار کاربر قرار می‌گیرد، اما در محاسبات نمادین با امکاناتی که درون متلب وجود دارد، مقدار دقیق با استفاده از فرمول‌ها بجای محاسبات بدست می‌آید.

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

دستور syms  که به دنبال آن نمادهای مورد نظر می‌آید، نمادهای مورد نظر کاربر را در متلب ثبت می‌کند. مثالی از این دستور:

≫ syms x,y,z

با وارد کردن این دستور، سه نماد x,y,z در فضای کار متلب، ثبت می‌شود و می‌توان از آن‌ها استفاده کرد، برای مثال با استفاده از نمادهای تعریف‌شده، تابع (sin(x را بصورت نمادین تعریف می‌کنیم.

≫ f = sin(x)

پس از اجرای دستور بالا، یک تابع نمادین با نام f تعریف می‌شود و مقدار آن برابر تابع نمادین (sin(x  خواهد بود. حال می‌توان با استفاده از دستورات و امکانات موجود برای توابع نمادین، مشتق، انتگرال و یا حد را برای این تابع و یا توابع تعریف‌شده‌ی دیگر محاسبه نمود.

 

محاسبه‌ی مشتق

برای محاسبه‌ی مقدار مشتق یک تابع در متلب، ابتدا لازم است تا متغیرهای نمادین مورد نیاز را تعریف نموده و سپس توابع را بصورت نمادین و با استفاده از نمادهای تعریف شده، تعریف کنیم. انجام محاسبه‌ی مشتق نیز بصورت نمادین است و نتیجه‌ای که متلب بعنوان خروجی در اختیار کاربر قرار می‌دهد نیز بصورت نمادین می‌باشد.

دستور diff مربوط به محاسبه‌ی مشتق (دیفرانسیل) یک تابع نمادین می‌باشد و پارامترهای ورودی آن به ترتیب نماد تابع تعریف‌شده، نماد متغیر مشتق و سپس مرتبه‌ی مشتق می‌باشد. کاربرد این دستور را در قالب یک مثال بررسی می‌کنیم:

≫ syms x

≫ f = sin(x)

≫ diff(f,x,1)

ans =

cos(x)

≫ f = x*exp(x)

≫ diff(f,x,5)

ans =

۵*exp(x) + x*exp(x)

 

محاسبه‌ی حد

یکی از کاربردهای محاسبات نمادین در متلب، محاسبه‌ی حد برای یک تابع نمادین است. برای این کار نیز ابتدا باید متغیرهای نمادین و توابع نمادین تعریف شود و سپس با استفاده از دستور limit مقدار حد را محاسبه نمود.

در دستور limit پارامترهای ورودی به ترتیب شامل نماد تابع مورد نظر، متغیر مورد حد در تابع و نقطه‌ی حدی می‌باشد. در مثال زیر به بررسی کاربرد این دستور می‌پردازیم:

≫ syms x

≫ f = (1 + 1/x) ^ x

≫ limit(f,x,inf)

ans =

exp(1)

 

≫ syms x

≫ f = (2*x^2 + 1) / (x^2 - 1)

≫ limit(f,x,inf)

ans =

۲

 

محاسبه‌ی انتگرال

محاسبه‌ی انتگرال در متلب نیز مانند حد و مشتق با استفاده از ویژگی محاسبات نمادین امکان‌پذیر خواهد بود. البته با توجه به اینکه انتگرال انواع روش‌ها و محدوده‌های متفاوتی دارد، دستورات متنوع‌تری برای محاسبه‌ی آن نیز وجود دارد.

برای بدست آوردن انتگرال نامعین یک تابع، با استفاده از دستور int که پارامترهای به ترتیب نماد تابع و نماد متغیر را دریافت می‌کند، به راحتی می‌توان حاصل را بدست‌آورد.

در انتگرال معین، دو قید برای متغیر انتگرال وجود دارد، لذا دستور انتگرال معین، همانند دستور انتگرال نامعین است، البته به اضافه‌ی دو پارامتر a بعنوان شروع محدوده و b بعنوان پایان محدوده‌ی انتگرال که به ترتیب پس از اعلام نماد متغیر وارد می‌شوند.

≫ syms x

≫ f = x*sin(x)

≫ int(f, x)"

ans =

sin(x) - x*cos(x)

≫ syms x

≫ f = exp(x) * (sin(x) + cos(x))

≫ int(f.x)

ans =

-exp(-x) * cos(x)

≫ syms x

≫ f = sin(sqrt(x))

≫ int(f, x, -2*pi, 2*pi)

ans =

?!

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

 

محاسبه‌ی مجموع سری

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

برای محاسبه‌ی مجموع یک سری، ابتدا متغیر نمادین و جمله‌ی عمومی نمادین آن را تعریف نموده و دستور symsum را با پارامترهای به ترتیب نماد جمله عمومی، نماد متغیر، اندیس شروع و اندیس پایان فراخوانی می‌کنیم.

همینطور برای محسابه‌ی سری تیلور یک سری، از دستور taylor و پارامترهای نماد سری، جمله‌ی اولیه و تعداد جمع‌ها استفاده می کنیم. در صورت ذکر نکردن جمله‌ی اولیه، سری مک‌لوران محاسبه می‌شود.

مثالی از مجموع سری و مجموع تیلور را مشاهده می‌کنیم:

≫ syms n

≫ f = n^2-n

≫ symsum(f, n, 1, inf);

≫ taylor(f, 1, inf);

 

حل معادلات دیفرانسیل

در معادلات دیفرانسیل، دستوراتی وجود دارد که به بررسی آن‌ها می‌پردازیم.

ابتدا برای حل یک معادله‌ی دیفرانسیل با دستور dsolve لازم است نمایش معادله را بصورت رشته با فرمت زیر بعنوان اولین پارامتر فراخوانی کنیم:

'{D{1}{2}={3' که در آن {۱}=مرتبه، {۲}=نماد و {۳}=معادله است.

بعنوان دومین ( یا بیشتر) پارامتر این دستور، باید شرایط معادله را به ترتیب بعنوان پارامتر در فراخوانی دستور اضافه کنیم. در مثال زیر، معادله‌ی dy/dx = 1+y^2 با شرایط اولیه‌ی y(0)=1 را محاسبه می‌کنیم:

≫ dsolve('Dy=1+y^2', 'y(0)=1')

ans =

tan(t+1/4*pi)

معادله‌ی d^2y/d^2x = cos(2x)-y با شرایط اولیه‌ی y(0)=1 و dy(0)/dx = 0

≫ dsolve('D2y=cos(2*x)-y', 'y(0)=1', 'Dy(0)=0', 'x')

ans =

(۱/۲*sin(x) + 1/6*sin(3*x))*sin(x) + (1/6*cos(3*x) – ۱/۲)

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

≫ simplify(ans)

ans =

-۲/۳*cos(x)^2 + 4/3*cos(x) + 1/3

 

تبدیل به لاپلاس و برعکس

برای محاسبه‌ی لاپلاس یک تابع یا معادله‌ی نمادین می‌توان از دستور laplace استفاده کرد. ابتدا متغیر نمادین و تابع نمادین را تعریف نموده و سپس نماد تابع یا متغیر را بعنوان پارامتر با دستور فراخوانی می‌کنیم. مثالی از این دستور:

≫ syms t

≫ f = exp(t)*cos(t)

≫ laplace(t)

ans =

(s-1)/((s-1)^2+1)

برای بدست آوردن معکوس لاپلاس یک نماد، دقیقاً مشابه قبل از دستور ilaplace استفاده می‌کنیم.

 

انتگرال‌گیری (غیر نمادی – عددی)

این دستورها منسوخ هستند و به زودی از متلب حذف خواهد شد و صرفاً جهت آشنایی توضیح داده می‌شود، برای کاربرد از دستورات جدید استفاده کنید.

برای انتگرال‌گیری یگانه بصورت عددی دو دستور وجود دارد که تنها در دقت محاسبه اختلاف دارند. دستور اول، دستور quad می‌باشد که پارامترهای آن به ترتیب نام تابع، حد بالا و حد پایین انتگرال است، در این دستور تابع باید بصورت “function” تعریف شود و سپس به شکل زیر مورد استفاده قرار گیرد:

| function y=f(x)

| y=sqrt((cos(x).^2+(sin(2.*x)).^2+x);

≫ quad(@f, 0, 2*pi);

دستور دیگر دقیقاً مانند دستور قبل عمل می‌کند با این تفاوت که یک پارامتر اضافه بعنوان دقت محاسبه در هنگام فراخوانی دریافت می‌کند. (مقدار پیش‌فرض آن ۱۰^-۶ است.) مثالی از دستور quadl را مشاهده می‌کنیم:

≫ quadl(@f, 0, 2*pi, 10^-8);

برای انتگرال‌گیری دوگانه از دستور dblquad استفاده می‌کنیم. در این دستور، پارامترها مانند دستور quadl می‌باشد، یعنی بصورت تقریبی انتگرال را محاسبه می‌کند، با این تفاوت که حد بالا و پایین را به ترتیب برای انتگرال داخلی و بیرونی دریافت می‌کند. مثالی از این دستور:

| function z=g(x.y)

| z=cos(x.*y)+x.*sin(y);

≫ dblquad(@g, -pi, 3*pi, 0, pi);

به همین صورت با دستور tripquad می‌توان انتگرال‌گیری سه‌گانه انجام داد.

| function w=h(x.y.z)

| w=(exp(x).*z.*sin(y)+x.*y.*z);

≫ tripquad(@h, 0, 1, -1, 2, 0, 1);

برای محاسبه‌ی انتگرال بصورت عددی، دستور integral در متلب وجود دارد. این دستور، انتگرال یگانه را محاسبه می‌کند و پارامترهای ورودی آن مانند دستور قدیمی quad است، با این تفاوت که برای تعریف تابع مورد نظر، بصورت inline عمل می‌کنیم. نگارش تعریف تابع بصورت inline را در مثال زیر می‌بینیم:

>> func = @(x) cos(sqrt(x));

>> integral(func, -pi/2, pi/2);

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

برای انتگرال‌های دوگانه و سه‌گانه نیز مانند دستورات dblquad و tripquad، دستورات integral2 و integral3 وجود دارد که کاربرد آن‌ها مانند integral است، با این تفاوت که پارامترهای بعدی آن، بازه‌های دوم و سوم را مشخص می‌کند. مثالی از کاربرد انتگرال دوگانه را با بازه‌ی متغیر بررسی می‌کنیم:

>> func = @(x,y) 1 ./ (sqrt(x+y)(1+x+y))

>> ymax = @(x) 1-x

>> integral2(func, 0, 1, 0, ymax)

 

معادلات دیفرانسیل (غیر نمادی – عددی)

متلب قادر است معادلات دیفرانسیل معمولی مقدار اولیه را حل کند. فرم کلی معادلات باید بصورت زیر باشد تا با تغییر متغیر، یک معادله‌ی مرتبه‌ی n را به n معادله‌ی مرتبه اول تبدیل کند.

f1

با SOLVER که پارامترهای آن به ترتیب نام تابع، بازه‌ی مورد نظر و مقادیر اولیه‌ی معادله است، می‌توان معادلات دیفرانسیل را حل نمود و خروجی این دستور یک ماتریس دو ستونه است که می‌توان به فرمی که در مثال آمده‌است، مقادیر x و y را مستقیماً دریافت کرد.

البته SOLVER مجموعه راه حل‌هایی است که متلب در اختیار کاربر قرار می‌دهد، تمام دستورات این مجموعه شامل موارد زیر است که هرکدام روش خاص خود را برای حل معادله دارد:

ode45, ode23, ode113, ode15s, ode23s, ode23t, ode23tb

برای مثال، معادله‌ی y''' - 3yy'' + y'sinx = 0 که در آن ۰ ≤ x ≤ ۱ و شرایط اولیه آن y(0)=0, y'(0)=-1, y''(0)=1 می‌باشد را با دستور ode45 حل کنیم.

ابتدا معادله مرتبه ۳ بالا را به سه معادله‌ی مرتبه‌ی اول تبدیل می‌کنیم:

y1' = y2

y2' = y3

y3' = 3*(y1)*(y3) - (y2)*sin(x)

y1(0) = 0

y2(0) = -1

y3(0) = 1

| function dy=f(x.y)

| dy = [y(2); y(3); 3*y(1)*y(3)-y(2)*sin(x)];

≫ [X Y] = ode45(@f, [0 1], [0; -1; 1]);

ستون اول ماتریس Y همان جواب معادله است.

در صورتی که بخواهیم جواب بدست آمده را رسم کنیم باید از دستور ezplot استفاده کنیم:

≫ ezplot(Y)

ode_plot_result

بسته به نوع معادلات که اصطلاحا به آنها سخت stiff و غیرسخت stiffness گفته می‌شود، روش حل آنها در MATLAB کمی متفاوت خواهد بود. اصطلاح سخت stiff برای آن دسته از معادلاتی بکار می‌رود که برای مثال در مقابل متغیر مستقلی همچون t چند متغیر وابسته مانند x و y و… وجود دارد، بگونه‌ای که اندازه‌ی مشتقات متغیرهای وابسته نسبت به متغیر مستقل، بطور قابل ملاحظه‌ای متفاوت است. در غیر اینصورت معادله غیرسخت نامیده می‌شود. همچنین معادلات سخت شامل آن دسته از معادلات دیفرانسیلی می‌شوند که حل آن‌ها با روش محاسبات عددی پایدار و همگرا نبوده و تنها راه حل آن‌ها، بسیار کوچک کردن گام در روش عددی می‌باشد. اگر در معادله دیفرانسیلی، متغیری وجود دارد که باعث تغییرات بسیار زیاد در جواب مساله می‌شود، این دسته را جز معادلات غیرسخت طبقه بندی می‌کنند.

کاربرد

دقت

نوع مساله

دستور

اکثر موارد (سعی شود جهت حل معادله ابتدا از این دستور استفاده شود)

متوسط

غیرسخت

ode45

حل مسائل دارای خطای خام (crude error)، حل مسائل تقریبا سخت

پایین

غیرسخت

ode23

حل مساله دارای خطای دقیق، معادلات مربوط به محاسبات عددی زمان بر

پایین تا بالا

غیرسخت

ode113

هنگامی که حل معادله با دستور ode45 بسیار کند پیش رود.

پایین تا متوسط

سخت

ode15s

مسائل دارای خطای خام (crude error) ، هنگامی که ماتریس جرم (mass matrix) ثابت باشد

پایین

سخت

ode23s

حل معادلات بدون میرایی عددی (numerical damping)

پایین

سخت

ode23t

حل مساله با خطای خام و معادله سخت

پایین

سخت

ode23tb

مثالی دیگر از معادلات دیفرانسیل:

| function dydt=vdp1(t,y)

| epsilon=5;

| w=2.466;

| f=5;

| dydt=[y(2) ; epsilon*(1-y(1)^2)*y(2)-y(1)+f*cos(w*t)];

>> [t,y]=ode45(@vdp1,[0 100],[1.2 0]);

>> xlabel('y1')

>> ylabel('y2')

>> plot(y(:,1),y(:,2))