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

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

Метод Адамса

Файл adams.f

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

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

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

Файл progon.f

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

Реализация на 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-го порядка