Лабораторна робота
Тема:
Розв’язування
систем
диференціальних рівнянь
Мета роботи: Вивчення методів розв’язування систем диференціальних рівнянь та
засобiв мови Сi для програмування процедур, оволодiння прийомами складання алгоритмiв i програм з пiдпрограмами.
Завдання:
вид дифрівняння | a3y'''+a2y''+a0y=b0
початкові умови | y=y'=y''=0 при t=0
b0 | C0k1k2
a0 | 1+b0
a1 | T1+T3
a2 | T2+T1T3
a3 | T1T2
C0k | a1a2/(a3k1k2)-1/(k1k2)
k1 | 3
k2 | 5
T1 | 0,1
T2 | 0,01
T3 | 0,2
tвст | 2,0
Короткі теоретичні відомості:
Розв’язування диференціального рівняння
Постановка задачі (задача Коші) має вигляд диференціального рівняння з початковими умовами
yР = f(t,y), y = y0 при t = t0. t Є [a, b].
Для її наближеного розв’язання застосовуються так звані однокрокові методи: Ейлера, Ейлера покращений, Ейлера-Коші та Рунге-Кута. Їх суть полягає в тому, що діапазон інтегрування [a, b] ділять на n елементарних відрізків довжиною h. Значення шуканої функції в точці t0=a відомо з початкових умов, а її обчислення в першій і наступних точках аж до точки tn=b виконують за поданими нижче формулами. При цьому h=(b-a)/n, t0 = a, tn = b, ti+1 = ti+h, yi=f(ti), i=0,1,2, ... n.
Алгоритм розв’язування дифрівняння однокроковими методами представляє собою цикл з накопиченням суми.
Метод Ейлера найпростіший, але він має лише перший порядок точності. Геометрично метод представляє собою перехід у кожну наступну точку, починаючи з заданої початковими умовами, вздовж дотичної до точного розв’язку, проведеної з попередньої точки. Цей метод враховує кут нахилу дотичної до кривої точного розв’язку лише в одній крайній точці елементарного відрізка. Ілюстрація першого кроку переходу в точку t1 подана на рисунку 2, де y1* – точне значення шуканої функції в точці t1, e – похибка. Формула Ейлера має такий вигляд: yi+1=yi+hf(xi,yi).
Метод Ейлера покращений має другий порядок точності. На кожному кроці обчислювального процесу знаходять проміжні значення шуканої функції yiп посередині кожного i-го елементарного відрізка (на відстані пів-кроку, тобто при ti+h/2) за методом Ейлера, а проміжні значення використовують для переходу в наступну точку. Геометрично це означає перехід у кожну наступну точку вздовж дотичної, проведеної до кривої точного розв’язку в середній точці. Формули мають такий вигляд:
yiп = yi + hf(ti, yi)/2;
yi+1 = yi + hf(ti+h/2, yiп).
Метод Ейлера-Коші теж має другий порядок точності. Для переходу в кожну наступну точку враховується кут нахилу дотичної до кривої точного розв’язку в обох крайніх точках (тобто в точках ti i ti+1) кожного елементарного відрізка. Для переходу в кожну наступну точку використовують середнє арифметичне значення тангенсів кутів нахилу дотичних до кривої точного розв’язку (похідних) в крайніх точках. Тут теж обчислюють проміжні значення. Формули мають такий вигляд:
yiп = yi + hf(ti, yi);
yi+1 = yi + h[f(ti+h, yiп) + f(ti, yi)]/2.
Метод Рунге-Кута має четвертий порядок точності. На кожному кроці циклічного обчислювального процесу знаходять допоміжні коефіцієнти a, b, c, d. Формули мають такий вигляд:
a = f(ti, yi);
b = f(ti+h/2, yi + ha/2);
c = f(ti+h/2, yi + hb/2);
d = f(ti+h, yi + hc);
yi+1 = yi + h(a + 2b + 2c + d)/6.
Перетворення заданого дифрівняння в систему дифрівнянь:
Для перетворення заданого дифрiвняння в систему рiвнянь першого по---рядку застосуємо таку пiдстановку:
Ваpiант 0: y=x3;
y'=x3'=x2;
y''=x3''=x2'=x1.
Система дифрiвнянь прийме вигляд
x1'=(b0-a0*x3-a1*x2-a2*x1)/a3;
x2'=x1;
x3'=x2.
Початковi умови: x1=x2=x3=0 при t=0. Вихiдна величина y=x3.
Таблиця ідентифікаторів
Ім’я
змінної | Ідентифікатор у програмі мовою
Сі |
Значення
y | yb[101] | Масив значень величини у (для графіка)
y[3] | Масив правих частин системи дифрівнянь
z[3] | Допоміжний масив правих частин системи дифрівнянь
f[3] | Масив правих частин системи дифрівнянь у підпрограмі
k1 | k1 | =3 – заданий коефіцієнт дифрівняння
k2 | k2 | =5 – заданий коефіцієнт дифрівняння
а0 | a0 | = 1+b0
а1 | a1 | =T1 + T3
а2 | a2 | =T2+T1T3
а3 | a3 | =T1T2
b0 | b0 | =C0k1k2
h | h | h=tвст/n – крок інтегрування
Cок | c0k | Коливна межа стійкості рішення
С1 | c1 | =0.2 – заданий коефіцієнт дифрівняння
tвст | twst | =2.0 час припинення коливань системи
Т1 | t1 | =0.1 – заданий коефіцієнт дифрівняння
Т2 | t2 | =0.01 – заданий коефіцієнт дифрівняння
Т3 | t3 | =0.2 – заданий коефіцієнт дифрівняння
i,j | Цілі змінні, параметри циклів
m | =3 – кількість дифрівнянь системи
n | =100 – число відрізків при інтегруванні
a[3],b[3],c[3], d[3] | Коефіцієнти формули Руннге-Кута
Програма:
#include<stdio.h>
float yb[101],y[3],z[3],f[3],a[3],b[3],c[3],d[3];
float k1=3,k2=5,a0,a1,a2,a3,b0,b1,h,c0;
float twst=2,t1=0.1,t2=0.01,t3=0.2;
int i,j,m=3,n=100,q;
Ejler_Pok(),Ejler_Kos(),Run_Kut(),systema();
/* Головна програма */
main()
{
h=twst/n;
a1=t1+t3;a2=t2+t1*t3;a3=t1*t2;
c0=a1*a2/(k1*k2*a3)-1/(k1*k2);
b0=c0*k1*k2;
a0=1+b0;
for(q=0;q<3;q++)
{
clrscr();
fprintf(stdout,"Результати розрахункiв:\n");
fprintf(stdout,"Коливне c0=%f\n",c0);
y[0]=y[1]=y[2]=0;yb[0]=0;
if(q==0){Ejler_Pok();printf("Метод Ейлера(покращений):\n");}
if(q==1){Ejler_Kos();printf("Метод Ейлера-Кошi:\n");}
if(q==2){Run_Kut(); printf("Метод Рунге-Кута:\n");}
for(i=1;i<=n;i+=5)
fprintf(stdout,"%f\t%f\t%f\t%f\t%f\n",
yb[i],yb[i+1],yb[i+2],yb[i+3],yb[i+4]);
bioskey(0);
}
}
/* Система дифрiвнянь */
systema()
{
f[0]=(b0-a0*y[2]-a1*y[1]-a2*y[0])/a3;
f[1]=y[0];
f[2]=y[1];
return;
}
/* Метод Ейлера покращений */
Ejler_Pok()
{
for(i=1;i<=n;i++)
{
for(j=0;j<m;j++)z[j]=y[j];
systema();
for(j=0;j<m;j++)y[j]=z[j]+h*f[j]/2;
systema();
for(j=0;j<m;j++)y[j]=z[j]+h*f[j];
yb[i]=y[2];
}
return;
}
/* Метод Ейлера-Кошi */
Ejler_Kos()
{
float a[4];
for(i=1;i<=n;i++)
{
for(j=0;j<m;j++)z[j]=y[j];
systema();
for(j=0;j<m;j++){a[j]=f[j];y[j]=z[j]+h*a[j];}
systema();
for(j=0;j<m;j++){y[j]=z[j]+h*(a[j]+f[j])/2;}
yb[i]=y[2];
}
return;
}
/* Метод Рунге-Кута */
Run_Kut()
{
for(i=1;i<=n;i++)
{
for(j=0;j<m;j++)z[j]=y[j];
systema();
for(j=0;j<m;j++){a[j]=f[j]*h;y[j]=z[j]+a[j]/2;}
systema();
for(j=0;j<m;j++){b[j]=f[j]*h;y[j]=z[j]+b[j]/2;}
systema();
for(j=0;j<m;j++){c[j]=f[j]*h;y[j]=z[j]+c[j];}
systema();
for(j=0;j<m;j++){d[j]=f[j]*h;y[j]=z[j]+(a[j]+2*(b[j]+c[j])+d[j])/6;}
yb[i]=y[2];
}
return;}
Результати:
Результати розрахункiв:
Коливне c0=0.533333
Метод Ейлера(покращений):
0.000000 0.054400 0.179968 0.374952 0.623338
0.899163 1.171071 1.406969 1.578456 1.664620
1.654802 1.550043 1.363049 1.116696 0.841252
0.570677 0.338444 0.173401 0.096192 0.116680
0.232685 0.430213 0.685144 0.966175 1.238665
1.468904 1.628282 1.696844 1.665790 1.538599
1.330646 1.067343 0.781036 0.507028 0.279232
0.125969 0.066452 0.108373 0.246917 0.465308
0.736839 1.028141 1.303281 1.528192 1.674894
1.724975 1.671896 1.521846 1.293022 1.013424
0.717426 0.441556 0.219983 0.080289 0.040035
0.104555 0.266250 0.505486 0.792971 1.093347
1.369533 1.587306 1.719538 1.749578 1.673346
1.499879 1.250270 0.955096 0.650676 0.374603
0.161102 0.036786 0.017344 0.105554 0.290902
0.550828 0.853475 1.161588 1.437105 1.645854
1.761794 1.770260 1.669819 1.472494 1.202325
0.892442 0.581005 0.306501 0.102992 -0.004113
-0.001228 0.111685 0.321068 0.601387 0.918250
1.232625 1.505646 1.703418 1.801227 1.786623 | Результати розрахункiв:
Коливне c0=0.533333
Метод Ейлера-Кошi:
0.000000 0.054400 0.179968 0.374952 0.623338
0.899163 1.171071 1.406969 1.578456 1.664620
1.654802 1.550043 1.363049 1.116696 0.841252
0.570677 0.338444 0.173401 0.096192 0.116680
0.232685 0.430213 0.685144 0.966175 1.238665
1.468904 1.628282 1.696844 1.665790 1.538599
1.330646 1.067343 0.781036 0.507028 0.279232
0.125969 0.066452 0.108373 0.246917 0.465308
0.736839 1.028141 1.303281
