Курсовая работа: 3 задание
Описание
Характеристики курсовой работы
Список файлов
- 3 задание
- zad3.txt 3,85 Kb
program bl1
use dflib
integer(4) color, paint
real(8) :: xmin=0, xmax=2, ymin=0, ymax=2, niu=0.1, theta=0.6, u0=1, eps=0.1**10, dt=0.01, f2=0.333984375
integer, parameter :: Nx=400, Ny=400
integer n, i, j, it
real(8) x(Nx), y(Ny), u(Nx,Ny), v(Nx,Ny), alpha(Ny), beta(Ny), A, B, C, D, ue(Nx,Ny), ve(Nx,Ny), ui, vi, du(Ny), delta, up(Ny), dx, dy, zeta(Nx,Ny), f(3), k1(3), k2(3), k3(3), k4(3), t, H(Ny),delta1
type(wxycoord) s
paint=setwindow(.true.,ymin,dble(0),u0,y max)
call init
call solve
contains
subroutine init
dx=(xmax-xmin)/(Nx-1.0); dy=(ymax-ymin)/(Ny-1.0)
do n=1, Nx
x(n)=dx*n
end do
do i=1, Ny
y(i)=dy*(i-1)
end do
do n=2, Nx
! точное решение
ue(1,:)=u0
do i=1, Ny
zeta(n,i)=y(i)*sqrt(u0/(x(n)*niu))
!if (zeta(n,i).gt.10) then
! zeta(n,i)=10
!end if
f(1)=0; f(2)=0; f(3)=f2
do t=0,zeta(n,i),dt
k1(1)=df1(t,f); k1(2)=df2(t,f); k1(3)=df3(t,f)
k2(1)=df1(t+0.5*dt,f+0.5*k1*dt); k2(2)=df2(t+0.5*dt,f+0.5*k1*dt); k2(3)=df3(t+0.5*dt,f+0.5*k1*dt)
k3(1)=df1(t+0.5*dt,f+0.5*k2*dt); k3(2)=df2(t+0.5*dt,f+0.5*k2*dt); k3(3)=df3(t+0.5*dt,f+0.5*k2*dt)
k4(1)=df1(t+dt,f+k3*dt); k4(2)=df2(t+dt,f+k3*dt); k4(3)=df3(t+dt,f+k3*dt)
f=f+dt/6*(k1+2*k2+2*k3+k4)
end do
ue(n,i)=u0*f(2); ve(n,i)=sqrt(niu*u0/x(n))*0.5*(zeta(n,i) *f(2)-f(1))
end do
end do
u(2:Nx,1)=0; v(:,:)=0; u(:,Ny)=u0; u(1,:)=u0
!H=ue(Nx,:)
!u=ue; v=ve
alpha(2)=0;
end subroutine init
subroutine solve
!схема
do n=1,Nx-1
delta=10; !write(*,*) n
it=0;
beta(2)=u(n+1,1)
do while (delta.gt.eps)
up(:)=u(n+1,:); delta1=0;
do i=2,Ny-1 ! прямая прогонка
if (it.gt.0) then
ui=u(n+1,i)*theta+u(n,i)*(1-theta)
else
ui=u(n,i)
end if
vi=v(n,i)
A=-theta*(0.5*vi*dx/dy+niu*dx/dy**2)
B=theta*(0.5*vi*dx/dy-niu*dx/dy**2)
C=-ui-theta*2*niu*dx/dy**2
D=-ui*u(n,i)+(1-theta)*(0.5*vi*(u(n,i+1) -u(n,i-1))*dx/dy-niu*(u(n,i+1)-2*u(n,i)+ u(n,i-1))*dx/dy**2)
alpha(i+1)=B/(C-alpha(i)*A); beta(i+1)=(A*beta(i)+D)/(C-alpha(i)*A)
end do
do i=Ny-1,1,-1 ! обратная прогонка
u(n+1,i)=u(n+1,i+1)*alpha(i+1)+beta(i+1)
end do
do i=1,Ny-1
du(i)=(u(n+1,i)-u(n,i))/dx
end do
do i=1,Ny
v(n,i)=0
do j=1,i-1 ! интеграл
v(n,i)=v(n,i)-0.5*dy*(du(j)+du(j+1))
end do
delta1=delta1+abs(u(n+1,i)-up(i));
end do
delta=delta1; !write(*,*) delta
it=it+1
end do
!write(*,*) it
end do
!write(*,*) u(:,:)
!call drawres(u,10,Ny+41)
!call drawres(v,Nx+20,Ny+41)
!call drawres(ue,10,2*Ny+50)
!call drawres(ve,Nx+20,2*Ny+50)
color=RGBToInteger(255,0,0)
call res(ue(200,:),y,color)
color=RGBToInteger(0,255,255)
call res(u(200,:),y,color)
end subroutine solve
subroutine drawres(T,xc,yc)
real(8) T(Nx,Ny),Tmin,Tmax,r0(2)
integer(4) color(1276)
integer xc, yc
r0(1)=xc
r0(2)=yc
!определение цвета
do i=1,1276
if (i.le.255) then
color(i)=RGBToInteger(256-i,0,255)
else if ((i.gt.255).and.(i.le.511)) then
color(i)=RGBToInteger(0,i-256,255)
else if ((i.gt.511).and.(i.le.766)) then
color(i)=RGBToInteger(0,255,255+511-i)
else if ((i.gt.766).and.(i.le.1021)) then
color(i)=RGBToInteger(i-766,255,0)
else if ((i.gt.1021).and.(i.le.1276)) then
color(i)=RGBToInteger(255,255-(i-1021),0 )
end if
end do
Tmax=maxval(T)
Tmin=minval(T)
!рисование решения
do i=1,Nx
do j=1,Ny
paint=SetPixelRGB(r0(1)+i,r0(2)-j,color( int(((T(i,j)-Tmin)/(Tmax-Tmin))*1275+1.5 )))
end do
end do
end subroutine drawres
subroutine res(x,y,color)
Integer(4) color
real(8) x(Ny),y(Ny)
paint=SetColorRGB(color)
call MoveTo_w(x(1),y(1),s)
do i=2,Ny
paint=LineTo_w(x(i),y(i))
end do
end subroutine res
real(8) function df1(x,f)
real(8) x,f(3)
df1=f(2)
end function df1
real(8) function df2(x,f)
real(8) x,f(3)
df2=f(3)
end function df2
real(8) function df3(x,f)
real(8) x,f(3)
df3=-0.5*f(1)*f(3)
end function df3
end
Старожил
