Автор работы: Пользователь скрыл имя, 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-го порядка
Метод Адамса
Метод прогонки
Сводная таблица результатов
Вывод
Список используемой литературы
Результат работы программы:
Файл 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- >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( >,yh2(i-3))-9*f(x-3*h,yh2(i-4) 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))/ return end |
Результат работы программы:
Файл 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;
Результат работы программы:
Файл 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;
Результат работы программы:
Файл 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(
x=x+h;
i=i+1;
end;
printf("Result:\ny(%1.1f)=%f\
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,
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))/
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(
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(
c(i)=1/(m(i)-n(i)*c(i-1));
d(i)=funF(x)*h*h-n(i)*c(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;
Результат работы программы:
Файл 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.
Результат работы программы:
Файл 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.
Результат работы программы:
Файл 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)/
k4:=h*Func(x+h/2,y+3*k1/8+9*
k5:=h*Func(x+h,y+3*k1/2-9*k3/
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])-
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]
x := x+h;
i := i + 1;
end;
writeLn('eps = ', ((yh[n]-yh2[2*n])/(exp(4*ln(2)
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))
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)+
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)+
c[i]:=1/(m[i]-n[i]*c[i-1]);
d[i]:=FunD(x)*h*h-n[i]*c[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.
Результат работы программы:
Файл 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
Результат работы программы:
Файл 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
Результат работы программы:
Информация о работе Методы решения обыкновенных дифференциальных уравнений