Методы решения обыкновенных дифференциальных уравнений

Автор работы: Пользователь скрыл имя, 01 Апреля 2012 в 19:49, лабораторная работа

Краткое описание

Постановка задачи 3
Ручной счёт и алгоритмы 3
Метод Эйлера 3
Метод Рунге-Кутты 4-го порядка 4
Метод Рунге-Кутты 5-го порядка 5
Метод Адамса 6
Правило Рунге для оценки погрешности 6
Метод прогонки 7
Реализация на языках программирования 9
Реализация на С++ 9
Метод Эйлера 9
Метод Рунге-Кутты 4-го порядка 10
Метод Рунге-Кутты 5-го порядка 11
Метод Адамса 12
Метод прогонки 14
Реализация на Fortran 15
Метод Эйлера 15
Метод Рунге-Кутты 4-го порядка 16
Метод Рунге-Кутты 5-го порядка 17
Метод Адамса 18
Метод прогонки 20
Реализация на SciLab 21
Метод Эйлера 21
Метод Рунге-Кутты 4-го порядка 22
Метод Рунге-Кутты 5-го порядка 23
Метод Адамса 23
Метод прогонки 25
Реализация на Pascal 26
Метод Эйлера 26
Метод Рунге-Кутты 4-го порядка 27
Метод Рунге-Кутты 5-го порядка 28
Метод Адамса 29
Метод прогонки 30
Реализация на Basic 32
Метод Эйлера 32
Метод Рунге-Кутты 4-го порядка 33
Метод Рунге-Кутты 5-го порядка 34
Метод Адамса 35
Метод прогонки 36
Сводная таблица результатов 38

Содержание работы

Постановка задачи
Ручной счёт и алгоритмы
Метод Эйлера
Метод Рунге-Кутты 4-го порядка
Метод Рунге-Кутты 5-го порядка
Метод Адамса
Правило Рунге для оценки погрешности
Метод прогонки
Реализация на языках программирования
Реализация на С++
Метод Эйлера
Метод Рунге-Кутты 4-го порядка
Метод Рунге-Кутты 5-го порядка
Метод Адамса
Метод прогонки
Реализация на Fortran
Метод Эйлера
Метод Рунге-Кутты 4-го порядка
Метод Рунге-Кутты 5-го порядка
Метод Адамса
Метод прогонки
Реализация на SciLab
Метод Эйлера
Метод Рунге-Кутты 4-го порядка
Метод Рунге-Кутты 5-го порядка
Метод Адамса
Метод прогонки
Реализация на Pascal
Метод Эйлера
Метод Рунге-Кутты 4-го порядка
Метод Рунге-Кутты 5-го порядка
Метод Адамса
Метод прогонки
Реализация на Basic
Метод Эйлера
Метод Рунге-Кутты 4-го порядка
Метод Рунге-Кутты 5-го порядка
Метод Адамса
Метод прогонки
Сводная таблица результатов
Вывод
Список используемой литературы

Содержимое работы - 1 файл

отчёт.doc

— 979.50 Кб (Скачать файл)

 

Результат работы программы:

Метод Адамса

Файл adams.f

1

2-5

6

9   72

 

 

 

program adams

f(x,y)= exp(x*y)

real::x0, xN, h, x, k1, k2, k3, k4, k5, y, y0

integer::i

integer::n

real, allocatable::yh(:)

real, allocatable::yh2(:)

print *,'Vvedite x0: '

read *,x0

print *,'Vvedite xN: '

read *,xN

print *,'Vvedite shag: '

read *,h

print *,'Vvedite y(x0): '

read *, y0

n = ceiling(abs(xN-x0)/h);

allocate(yh(0:n),yh2(0:2*n))

x = x0; y = y0

i = 0

do while(x<=xN)

              yh(i) = y

k1 = h*f(x,y)

k2 = h*f(x + h/2.0, y + k1/2.0)

k3 = h*f(x + h/2.0, y + k2/2.0)

k4 = h*f(x + h, y + k3)

y = y+(k1+2.0*k2+2.0*k3+k4)/6.0

              i = i+1

x = x+h

end do

i = 4

x = x0 + i*h

do while(x<=xN)

              yh(i)=yh(i-1)+h*(55*f(x,yh(i-1))-59*f(x-h,yh(i-2))+37*f(x-2*h,yh(

>i-3))-9*f(x-3*h,yh(i-4)))/24;

              x=x+h

              i=i+1

end do

print *, 'Result: '

write(*,'(a,f3.1,a,f8.6)') 'y(', xN,')=', yh(n)

 

y = y0; x = x0; h = h/2.0

i = 0

do while(x<=xN)

              yh2(i) = y

k1 = h*f(x,y)

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

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

k4 = h*f(x + h, y + k3)

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

              i = i+1

x = x+h

end do

i = 8

x = x0 + i*h

do while(x<=xN)

              yh2(i)=yh2(i-1)+h*(55*f(x,yh2(i-1))-59*f(x-h,yh2(i-2))+37*f(x-2*h

>,yh2(i-3))-9*f(x-3*h,yh2(i-4)))/24;

              x = x+h

              i = i+1

end do

 

write(*,'(a, f8.6)') 'eps = ', (yh(n)-yh2(2*n))/(2**4-1)

end program adams

 

Результат работы программы:

где eps – погрешность вычисленного решения.

Метод прогонки

Файл progon.f

1

2-5

6

9   72

 

 

 

program progonka

              real :: a, b, x, h, y(100) ,c(100), d(100), m(100), n(100)

              integer :: i,K

              print *, 'Enter a = '

read (*,*) a

              print *, 'Enter b = '

read (*,*) b

              if (a >= b) then

stop 'bad interval'

              end if

              print *, 'Enter N = '

read (*,*) K

              print *, 'Enter y(',a, ') = '

read (*,*) y(1)

              print *, 'Enter y(',b, ') = '

read (*,*) y(K+1)

              h = (b-a)/K

              m(1) = -2 + funB(a)*h/funA(a)

              n(1) = 1 - funB(a)*h/funA(a) + funC(a)*h*h/funA(a)

              c(1) = 1/m(1)

              d(1) = -n(1)*y(1) + funF(a)*h*h

              do i = 2,K+1

                x = a + (i-1)*h

                m(i) = -2 + funB(x)*h/funA(x)

                n(i) = 1 - funB(x)*h/funA(x) + funC(x)*h*h/funA(x)

                c(i) = 1/(m(i) - n(i)*c(i-1))

                d(i) = funF(x)*h*h - n(i)*c(i-1)*d(i-1)

end do

              do i = K+1,3,-1

                y(i-1) = c(i-2)*(d(i-2) - y(i))

end do

              print *,'              x                            y'

do i = 1, K+1

                print *, a+(i-1)*h,' ',y(i)

end do

  end program lab4_progonka

 

real function funA(x)

real :: x

funA = 0.4+2*x+x*x/2

return

end

 

real function funB(x)

real :: x

funB = x**0.7

return

end

 

real function funC(x)

real :: x

funC = -0.95*x*x

return

end

 

real function funF(x)

real :: x

funF = (sqrt(x)+cos(2*x/3.1416))/funA(x)

return

end

 

Результат работы программы:

 

Реализация на SciLab

Метод Эйлера

Файл eiler.sce

function[res]=func(x,y);

  res=exp(x*y);

endfunction;

disp('Method eilera');

disp('y''=exp(x*y)');

x0=input('Vvedite x0: ');

xN=input('Vvedite xN: ');

h=input('Vvedite shag: ');

y=input('Vvedite y0: ');

printf("  x   y  dy\n");

while(x0<xN+h)

  dy=h*func(x0,y);

  printf("%1.1f%10f%10f\n", x0,y,dy);

  y=y+dy;

  x0=x0+h;

end;

 

Результат работы программы:

 

Метод Рунге-Кутты 4-го порядка

 

Файл RK4.sce

function[res]=func(x,y);

  res=exp(x*y);

endfunction;

x0=input('Vvedite x0: ');

xN=input('Vvedite xN: ');

h=input('Vvedite shag: ');

y=input('Vvedite y0: ');

x=x0;

printf("x\ty\n");

while(x<xN+h)

  printf("%f\t%f\n", x,y);

  k1 = h*func(x,y);

  k2 = h*func(x + h/2, y + k1/2);

  k3 = h*func(x + h/2, y + k2/2);

  k4 = h*func(x + h, y + k3);

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

  x=x+h;

end;

 

Результат работы программы:

Метод Рунге-Кутты 5-го порядка

 

Файл RK5.sce

function[res]=func(x,y);

  res=exp(x*y);

endfunction;

x0=input('Vvedite x0: ');

xN=input('Vvedite xN: ');

h=input('Vvedite shag: ');

y=input('Vvedite y0: ');

x=x0;

printf("x\ty\n");

while(x<xN+h)

  printf("%f\t%f\n", x, y);

  k1 = h*func(x,y)/3;

  k2 = h*func(x + h/3, y + k1)/3;

  k3 = h*func(x + h/3, y + k1/2 + k2/2)/3;

  k4 = h*func(x + h/2, y + 3*k1/8 + 9*k3/8)/3;

  k5 = h*func(x + h, y + 3*k1/2 - 9*k3/2 + 6*k4)/3;

  y=y+(k1 + 4*k4 + k5)/2;

  x=x+h;

end;

 

Результат работы программы:

 

Метод Адамса

Файл adams.sce

function[res]=func(x,y);

  res=exp(x*y)

endfunction;

x0 = input('Vvedite x0: ');

xN = input('Vvedite xN: ');

h = input('Vvedite shag: ');

i = 1;

y0 = input('Vveidte y(x0): ');

n = ceil(abs(xN-x0)/h);

y = y0; x = x0;

while(x<=xN)

  yh(i) = y;

  k1 = h*func(x,y);

  k2 = h*func(x + h/2, y + k1/2);

  k3 = h*func(x + h/2, y + k2/2);

  k4 = h*func(x + h, y + k3);

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

  x=x+h;

  i = i + 1;

end;

i = 5;

x = x0 + (i-1)*h;

while(x<=xN)

  yh(i)=yh(i-1)+h*(55*func(x,yh(i-1))-59*func(x-h,yh(i-2))+37*func(x-2*h,yh(i-3))-9*func(x-3*h,yh(i-4)))/24;

  x=x+h;

  i=i+1;

end;

printf("Result:\ny(%1.1f)=%f\n", xN, yh(n+1));

y = y0; x = x0; h = h/2;

while(x<=xN)

  yh2(i) = y;

  k1 = h*func(x,y);

  k2 = h*func(x + h/2, y + k1/2);

  k3 = h*func(x + h/2, y + k2/2);

  k4 = h*func(x + h, y + k3);

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

  x=x+h;

  i = i + 1;

end;

i = 9;

x = x0 + (i-1)*h;

while(x<=xN)

  yh2(i)=yh2(i-1)+h*(55*func(x,yh2(i-1))-59*func(x-h,yh2(i-2))+37*func(x-2*h,yh2(i-3))-9*func(x-3*h,yh2(i-4)))/24;

  x=x+h;

  i=i+1;

end;

printf("eps = %f\n", (yh(n+1)-yh2(2*n+1))/(2^4-1));

 

Результат работы программы:

где eps – погрешность вычисленного решения.

 

Метод прогонки

Файл progon.sce

function [res]=funA(x)

  res=0.4+2*x+x*x/2;

endfunction;

function [res]=funB(x)

  res=x^0.7;

endfunction;

function [res]=funC(x)

  res=-0.95*x*x;

endfunction;

function [res]=funF(x)

  res=(sqrt(x)+cos(2*x/3.1416))/funA(x);

endfunction;

a=input('vvedite a');

b=input('vvedite b');

N=input('vvedite n');

y(1)=input('vvedite y(0)');

y(N+1)=input('vvedite y(n)');

h=(b-a)/N;

m(1)=-2+funB(a)*h/funA(a);

n(1)=1-funB(a)*h/funA(a)+funC(a)*h*h/funA(a);

c(1)=1/m(1);

d(1)=-n(1)*y(1)+funF(a)*h*h;

for i=2:N

              x=a+(i-1)*h;

              m(i)=-2+funB(x)*h/funA(x);

              n(i)=1-funB(x)*h/funA(x)+funC(x)*h*h/funA(x);

              c(i)=1/(m(i)-n(i)*c(i-1));

              d(i)=funF(x)*h*h-n(i)*c(i-1)*d(i-1);

end;

for i=N+1:-1:3

                y(i-1)=c(i-2)*(d(i-2)-y(i));

end;

printf("x\t\ty\n");

for i=1:N+1

  printf("%f  %f\n",a+(i-1)*h,y(i));

end;

 

Результат работы программы:

Реализация на Pascal

Метод Эйлера

Файл eiler.pas

var

x0,xN,h,y,dy: Real;

 

Function Func(x,y: Real): Real;

begin

  Func:= exp(x*y);

end;

 

begin

WriteLN('Method eilera');

WriteLN('y''=exp(x*y)');

Write('Enter x0: ');

ReadLn(x0);

Write('Enter xN: ');

ReadLn(xN);

Write('Enter shag: ');

ReadLn(h);

Write('Enter y(0): ');

ReadLn(y);

WriteLn('x ydy');

While x0<xN+h do begin

  dy:=h*Func(x0,y);

  WriteLn(x0:1:5,' ',y:1:5,' ',dy:1:5);

  y:=y+dy;

  x0:=x0+h;

end;

end.

 

Результат работы программы:

Метод Рунге-Кутты 4-го порядка

Файл RK4.pas

var

x,x0,xN,h,y,k1,k2,k3,k4: Real;

 

Function Func(x,y: Real): Real;

begin

  Func:= exp(x*y)

end;

 

begin

WriteLn('Enter x0:');

ReadLn(x0);

WriteLn('Enter xN:');

ReadLn(xN);

WriteLn('Enter step:');

ReadLn(h);

WriteLn('Enter y(0):');

ReadLn(y);

x:=x0;

WriteLn('x              y');

While x < xN+h do begin

  WriteLn(x:1:1,'              ', y:1:5);

  k1:=h*Func(x,y);

  k2:=h*Func(x+h/2,y+k1/2);

  k3:=h*Func(x+h/2,y+k2/2);

  k4:=h*Func(x+h,y+k3);

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

  x:=x+h;

end;

end.

 

Результат работы программы:

Метод Рунге-Кутты 5-го порядка

Файл RK5.pas

var

x,x0,xN,h,y,k1,k2,k3,k4,k5: Real;

 

Function Func(x,y: Real): Real;

begin

  Func:=exp(x*y);

end;

 

begin

WriteLn('Enter x0:');

ReadLn(x0);

WriteLn('Enter xN:');

ReadLn(xN);

WriteLn('Enter step:');

ReadLn(h);

WriteLn('Enter y(0):');

ReadLn(y);

x:=x0;

WriteLn('x              y');

While x<xN+h do begin

  WriteLn(x:1:1, '              ', y:1:5);

  k1:=h*Func(x,y)/3;

  k2:=h*Func(x+h/3,y+k1)/3;

  k3:=h*Func(x+h/3,y+k1/2+k2/2)/3;

  k4:=h*Func(x+h/2,y+3*k1/8+9*k3/8)/3;

  k5:=h*Func(x+h,y+3*k1/2-9*k3/2+6*k4)/3;

  y:=y+(k1+4*k4+k5)/2;

  x:=x+h;

end;

end.

 

Результат работы программы:

Метод Адамса

Файл adams.pas

var

x,x0,xN,h,y0, y,k1,k2,k3,k4: Real;

i:Integer;

n:Integer;

yh:array of Real;

yh2:array of Real;

Function Func(x,y: Real): Real;

begin

  Func := exp(x*y);

end;

begin

WriteLn('Enter x0:');

ReadLn(x0);

WriteLn('Enter xN:');

ReadLn(xN);

WriteLn('Enter step:');

ReadLn(h);

n:= Round(abs(xN-x0)/h);

SetLength(yh,n);

SetLength(yh2,2*n);

WriteLn('Enter y(x0):');

ReadLn(y0);

x := x0; y := y0; i:= 0;

While x <= xN do begin

  yh[i] := y;

  k1:= h*Func(x,y);

  k2:= h*Func(x+h/2,y+k1/2);

  k3:= h*Func(x+h/2,y+k2/2);

  k4:= h*Func(x+h,y+k3);

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

  x:= x+h;

  i := i+1;

end;

i := 4;

x:= x0+i*h;

While x <= xN do begin

  yh[i]:= yh[i-1]+h*(55*Func(x,yh[i-1])-59*Func(x-h,yh[i-2])+37*Func(x-2*h,yh[i-3])-9*Func(x-3*h,yh[i-4]))/24;

  x:= x+h;

  i := i + 1;

end;

WriteLn('Result: ');

WriteLn('y(', xN:1:1, ')=', yh[n]:8:6);

y := y0; x := x0; h := h/2;

i := 0;

While x <= xN do begin

  yh2[i] := y;

  k1 := h*Func(x,y);

  k2 := h*Func(x+h/2,y+k1/2);

  k3 := h*Func(x+h/2,y+k2/2);

  k4 := h*Func(x+h,y+k3);

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

  x := x+h;

  i := i+1;

end;

i := 8;

x := x0+i*h;

While x <= xN do begin

  yh2[i] := yh2[i-1]+h*(55*Func(x,yh2[i-1])-59*Func(x-h,yh2[i-2])+37*Func(x-2*h,yh2[i-3])-9*Func(x-3*h,yh2[i-4]))/24;

  x := x+h;

  i := i + 1;

end;

writeLn('eps = ', ((yh[n]-yh2[2*n])/(exp(4*ln(2))-1)):8:6);

end.

 

Результат работы программы:

где eps – погрешность вычисленного решения.

Метод прогонки

Файл progon.pas

var

a,b,x,h: Real;

N1,i:Integer;

y,c,d,m,n:array of Real;

 

Function FunA(x: Real): Real;

begin

  FunA:=0.4+2*x+x*x/2

end;

 

Function FunB(x: Real): Real;

begin

  FunB:=exp(0.7*ln(x));

end;

 

Function FunC(x: Real): Real;

begin

  FunC:=-0.95*x*x;

end;

 

Function FunD(x: Real): Real;

begin

  FunD:=(sqr(x)+cos(2*x/3.1416))/FunA(x);

end;

 

 

begin

WriteLn('Enter a:');

ReadLn(a);

WriteLn('Enter b:');

ReadLn(b);

WriteLn('Enter N:');

ReadLn(N1);

 

SetLength(y,N1+1);

SetLength(c,N1+1);

SetLength(d,N1+1);

SetLength(m,N1+1);

SetLength(n,N1+1);

 

WriteLn('Enter y(',a:1:5,')=');

ReadLn(y[0]);

WriteLn('Enter y(',b:1:5,')=');

ReadLn(y[N1]);

h:=(b-a)/N1;

 

m[0]:=-2+FunB(a)*h/FunA(a);

n[0]:=1-FunB(a)*h/FunA(a)+FunC(a)*h*h/FunA(a);

c[0]:=1/m[0];

d[0]:=-n[0]*y[0]+FunD(a)*h*h;

 

for i:=1 to N1 do begin

  x:=a+i*h;

  m[i]:=-2+FunB(x)*h/FunA(x);

  n[i]:=1-FunB(x)*h/FunA(x)+FunC(x)*h*h/FunA(x);

  c[i]:=1/(m[i]-n[i]*c[i-1]);

  d[i]:=FunD(x)*h*h-n[i]*c[i-1]*d[i-1];

end;

 

for i:=N1 downto 2 do

  y[i-1]:=c[i-2]*(d[i-2]-y[i]);

Writeln('x              y');

for i:=0 to N1 do

  WriteLn((a+i*h):1:5, '  ', y[i]:1:5);

 

end.

 

Результат работы программы:

Реализация на Basic

Метод Эйлера

Файл eiler.bas

function f(x as Double, y as Double) as Double

  f= exp(x*y)

end function

rem eiler

  dim x0 as Double,xN as Double

  dim h as Double,y as Double

  dim dy as Double

  dy = 0.0

  print "Method eilera"

  print "y'=exp(x*y)"

  input "Enter x0: ", x0

  input "Enter xN: ", xN

  input "Enter shag: ", h

  input "Enter y(0): ", y

  print "   x    y   dy"

  do while x0<xN+h

    dy=h*f(x0,y)

    print using "##.#  ##.######  ##.######"; x0; y; dy

    y = y + dy

    x0 = x0 + h

  loop

end

 

Результат работы программы:

 

Метод Рунге-Кутты 4-го порядка

Файл RK4.bas

function f(x as Double, y as Double) as Double

  f= exp(x*y)

end function

rem RK4

  dim x0 as Double,xN as Double

  dim h as Double

  dim y as Double

  dim k1 as Double, k2 as Double, k3 as Double, k4 as Double

  dim x as Double

  input "Enter x0: ", x0

  input "Enter xN: ", xN

  input "Enter step: ", h

  input "Enter y(0): ", y

  x = x0

  print "x              y"

  do while x < xN+h

    print using "#.#              ##.######"; x; y

    k1 = h*f(x,y)

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

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

    k4 = h*f(x + h, y + k3)

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

    x = x + h

  loop

end

Результат работы программы:

Метод Рунге-Кутты 5-го порядка

Информация о работе Методы решения обыкновенных дифференциальных уравнений