Розв’язування звичайних диференціальних рівнянь і систем методами Ойлера та Рунге-Кутта

Виконав: ст.гр.АВ-11 Скібель П.

Лабораторна  робота №13

«Розв’язування звичайних диференціальних рівнянь

і систем методами Ойлера та Рунге-Кутта»

Варіант №14

Завдання 1. Розв’язати звичайне диференціальне рівняння [pic 1] із заданою початковою умовою y(0)=0 на інтервалі x∈[0; 2] методами:  
1) Ойлера;  2) Рунге-Кутта; 3) в Matlab за допомогою функції ode45().

1. Метод Ойлера

Формула методу Ойлера для числового розв’язування ЗДР:

[pic 2].

[pic 3]

Програма в С

#include<stdio.h>

#include<math.h>

float f_lab13(float x, float y)  {return exp(-y)+cos(8*pow(x,2));}

int main()

            {int n=10;  float  xa=0, xb=2, h, x0, y0=0, y1;   h=(xb-xa)/n;

              for(x0=xa; x0<=xb+h/3; x0+=h)

                  {printf("\n x0=%.1f     y0=%f ", x0, y0);

                    y1=y0+h*f_lab13(x0,y0);

                    y0=y1;

                  }

              }

Таблиця 1

Результати виконання програми

Програма в Matlab

clc; clear all;

xa=0; xb=2; n=10;  h=(xb-xa)/n; y0=0;

for x0=xa:h:xb;

    disp('    x0        y0'); disp([x0,y0]);

    y1=y0+h*f_lab13(x0,y0);

    y0=y1;

end

Файл-функція:

function z=f_lab13(x,y)

         z=exp(-y)+cos(8*power(x,2));

Таблиця 2

Результати виконання програми

2. Метод Рунге-Кутта

Формули методу Рунге-Кутта для числового розв’язування ЗДР:

[pic 4]

[pic 5]

[pic 6]

Програма в С

#include<stdio.h>

#include<math.h>

float f_lab13(float x, float y)  {return  exp(-y)+cos(8*pow(x,2));}

int main()

            {int  n=10;  float  xa=0, xb=2, h, x0, y0=0, y1, k1, k2, k3, k4;

               h=(xb-xa)/n;

               for(x0=xa; x0<=xb+h/3; x0+=h)

                  {printf("\n x0=%.1f y0=%f",x0,y0);

                      k1=h*f_lab13(x0, y0);

                      k2=h*f_lab13(x0+h/2, y0+k1/2);

                      k3=h*f_lab13(x0+h/2, y0+k2/2);

                      k4=h*f_lab13(x0+h, y0+k3);

                      y1=y0+(k1+2*k2+2*k3+k4)/6;

                      y0=y1;

                   }

              }

Таблиця 3

Результати виконання програми:

Програма в Matlab

clc; clear all;

xa=0; xb=2; n=10;

h=(b-a)/n; y0=0;

for x0=xa:h:xb;

    disp('    x0        y0'); disp([x0,y0]);

    k1=h*f_lab13(x0,y0);

    k2=h*f_lab13(x0+h/2,y0+k1/2);

    k3=h*f_lab13(x0+h/2,y0+k2/2);

    k4=h*f_lab13(x0+h,y0+k3);

    y1=y0+(k1+2*k2+2*k3+k4)/6;