Metod konačnih razlika u Fortranu

Metod konačnih razlika u Fortranu 22.02.2010. u 10:45 · pre 203 meseca #1

Naime, imam problem sa programiranjem u Fortranu koji je sledeće prirode. Treba numerički rešiti parcijalnu diferencijalnu jednačinu metodom kočnanih razlika. Jednačina je prikačena.

Kod koji sam napisao je:

Code:
PROGRAM NUMERICKO

  dimension pr(-1:100,0:300000)

  REAL theta,DZ,D,A,c,indeks,Z,korakuglovni, ugao, kriticni
  INTEGER i,j,n,m
  DOUBLE PRECISION  w,Dtheta
  PARAMETER (c=299792458, w=314, korakuglovni=0.5, DZ=0.001, kriticni=20,D=0.000117, A=0,4, indeks=1.5 )

  PRINT *,'Unesite duzinu vlakna'
  READ (*,*)Z

  PRINT *,'Unesite ugao pod kojim se ubacuje signal u stepenima u stepenima' !ugao mora da bu edmanji od 15 stepeni
  READ (*,*)ugao

  Dtheta=korakuglovni*0.017453292
  theta=ugao*0.017453292
  thetaC=kriticni*0.017453292

  m=Z/DZ
  n=kriticni/korakuglovni
  PRINT *,n,m

  open(5,FILE='popuna1.dat')

  !unos impulsa na pocetku vlakna

  DO i = 0, n, 1
    DO j = 0, m, 1

         if (((i*korakuglovni).EQ.ugao .AND. (j.EQ.0))) then

            pr(i,j) = 1

        else

            pr(i,j) = 0

        end if
      END DO
   END DO

  open(6,FILE='racun1.dat')

  DO j = 0, m-1, 1
    DO i = -1, n+1, 1

     if ((i .LT. 0) .OR. (i.GT.n-1)) then

            pr(i,j)=0
    else

     if (i .EQ. 0) then

         p1=(2.*DZ)/(Dtheta**2.)
         p2=(4.*DZ)/(Dtheta**2.)-1.
         p3=(2.*DZ)/(Dtheta**2.)
            pr(i,j+1)=p1*pr(i-1,j)-p2*pr(i,j)+p3*pr(i+1,j)

     else

     if (i .EQ. 1) then

          pr(i,j)=pr(0,j)

         else

                p4=(DZ*D)/(Dtheta**2.)-(DZ*D)/(2.*Dtheta*i*Dtheta)
                p5=1.-(2.*DZ*D)/(Dtheta**2.)-DZ*A*((Dtheta*i)**2.)
                p6=(DZ*D)/(2.*Dtheta*i*Dtheta)+(DZ*D)/(Dtheta**2.)
                pr(i,j+1)=p4*pr(i-1,j)+p5*pr(i,j)+p6*pr(i+1,j)

            end if
        end if
      end if
      END DO
  END DO

  open(7,FILE='ispis1.dat')

   DO j = m-2, m, 1
      DO i = 0, n, 1
     if (j .EQ. m) then
         write (7,*)i,pr(i,j)
     end if
     END DO
  END DO

END PROGRAM

problem je u tome što funkcija treba da daje raspodelu koja se pri većim dužinama pomera ka 0 (na X osi), tj. maksimum treba da je na 0, a meni funkcija ima približni oblik Gausove raspodele oko vrednosti ugla pod kojim je ubačen signal. Znam da je problem negde u programiranju, ali ne znam gde. Prva tri uslova u glavnoj petlji pretstavljaju granične uslove.

Ako neko može da pomogne bio bih veoma zahvalan.

Ako je potreban još neki podatak tu sam

jednacina.jpg29.03k

Re: Metod konačnih razlika u Fortranu 23.02.2010. u 11:34 · pre 203 meseca #2

Trebalo bi da navedeš jednačinu i granične uslove, sa očekivanim opsegom ulaznih konstanti, i kakvu diskretizacionu šemu primenjuješ (trenutno si zakačio samo krajnju diskretizaciju).

Da ne bi skenirao rukopisne formule, vidi http://www.elitesecurity.org/t35291-Sve-La-TeX-na-ovom-forumu .

Re: Metod konačnih razlika u Fortranu 24.02.2010. u 10:52 · pre 203 meseca #3

Jednačina je: $\frac{\partial p(\theta,z)}{\partial z}=-A\theta^2p(\theta,z)+\frac{D}{\theta}\frac{\partial}{\partial\theta}(\theta\frac{\partial p(\theta,z)}{\partial\theta})$

što kad se sredi daje: $\frac{\partial p(\theta,z)}{\partial z}=-A\theta^2p(\theta,z)+\frac{D}{\theta}\frac{\partial p(\theta,z)}{\partial\theta}+D\frac{\partial^2p(\theta,z)}{\partial\theta^2}$

gde je: A-faktor multiplikacije drugog reda u izrazu za koeficijent gubitka snage
D-koeficijent sprezanja
z-dužina vlakna
θ-ugao prostiranja svetlosti u odnosu na osu vlakna

Granični uslovi su: $p(\theta_c,z)=0 \quad i \quad \frac{\partial p(\theta,z)}{\partial\theta}_{\theta=0} =0$

Primenom eksplicitnog metoda konačnih razlika, koristeći:

1) šemu centralne razlike za izvode $\frac{\partial p(\theta,z)}{\partial\theta}\quad i \quad \frac{\partial^2p(\theta,z)}{\partial\theta^2}$ dobijamo:

$(\frac{\partial p(\theta,z)}{\partial\theta})_{i,j}=\frac{p_{i+1,j}-p_{i-1,j}}{2\Delta\theta}$

$(\frac{\partial^2 p(\theta,z)}{\partial\theta^2})_{i,j})=\frac{p_{i+1,j}-2p_{i,j}+p_{i-1,j}}{\Delta\theta^2}$

2) šemu prednje razlike za izvod $\frac{\partial p(\theta,z)}{\partial z}$ dobijamo:

$(\frac{\partial p(\theta,z)}{\partial z})_{i,j}=\frac{p_{i,j+1}-p_{i,j}}{\Delta z}$

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

Re: Metod konačnih razlika u Fortranu 24.02.2010. u 11:49 · pre 203 meseca #4

dalje jednačina 2 postaje:

$p_{i,j+1}=(\frac{\Delta z D}{\Delta\theta^2}-\frac{\Delta z D}{2\Delta\theta\theta_{i,j}})p_{i-1,j} + (1-\frac{2\Delta z D}{\Delta\theta^2}-\Delta zA\theta_{i,j}^2)p_{i,j} + (\frac{\Delta z D}{2\Delta\theta\theta_{i,j}}+\frac{\Delta z D}{\Delta\theta^2})p_{i+1,j}$

pri čemu granični uslovi postaju: $p_{N,j}=0 \quad i \quad p_{0,j}=p_{1,j}$

Da bismo izbegli problem singularnosti u tačkama θ=0 koristimo: $\lim_{\theta \to 0} {\frac {1}{\theta}\frac{\partial}{\partial\theta}(\theta\frac{\partial p(\theta,z)}{\partial\theta}) = 2\frac{\partial^2p(\theta,z)}{\partial^2\theta} _{za\theta=0} }$

Upadni snop svetlosti je u obliku ravanskog talasa dat pomoću Dirakove delta funkcije: $p(\theta,0)=1\quad za \quad \theta=\theta_0\quad i \quad \cr p(\theta,0)=0\quad za \quad \theta\ne \theta_0$

gde je $\theta_0$ upadni ugao svetlosti u odnosu na osu vlakna.

Re: Metod konačnih razlika u Fortranu 24.02.2010. u 13:17 · pre 203 meseca #5

Trenutno me kopkaju dve stvari.

Prvo, u uslovu:

Code:
if (((i*korakuglovni) .EQ. ugao ...

pošto se radi o realnim vrednostima, ne vidim kako može da se pogodi tačno $i \Delta\theta = \theta_0$ . Ali, ako se ne pogađa, to bi značilo da program zapravo ne uspostavlja početni uslov, već da računa sa nekim „šumom“. Kako god bilo, ovaj uslov treba da bude:

Code:
abs(i * korakuglovni - ugao) .lt. Dtheta

(Takođe nema razloga da postavljaš $z$ -slojeve posle prvog na nulu, pošto je jednačina parabolička po $z$ .)

Drugo, nije mi jasna ova vratolomija sa levim graničnim uslovom. Koje je opravdanje za $\theta$ -sloj -1 fiksiran na nulu (pa onda računanje sloja 0 preko limita)? Ja bih rekao da ovde odgovara prost Nojmanov uslov, i da se petlje svode na:

Code:
! Početni uslov.
do i = 0, n
    if (abs(i * korakuglovni - ugao) .lt. Dtheta) then
        pr(i, 0) = 1.0
    else
        pr(i, 0) = 0.0
    end if
end do

! Marširanje po z.
do j = 0, m - 1
    pr(0, j) = pr(1, j)  ! levi granični uslov, Nojmanov
    pr(n, j) = 0.0       ! desni granični uslov, Dirihleov
    do i = 1, n - 1
        p4 = ...
        p5 = ...
        p6 = ...
        p(i, j + 1) = ...
    end do
end do

(I levi i desni granični uslov za tekući $z$ -sloj treba postaviti pre unutrašnje petlje.)

Ako i sa ove dve izmene ne radi kako je očekivano, onda treba ispitati uslove stabilnosti šeme. Pošto je jednačina nelinearna, ali „skoro“ linearna, onda bi trebalo da može Fon Nojmanova analiza za reprezentativne linearizacije (npr. pri $\theta = \Delta\theta$ i $\theta = \theta_c$ ) i date koeficijente $A$ i $D$ .

Re: Metod konačnih razlika u Fortranu 24.02.2010. u 13:46 · pre 203 meseca #6

Hvala puno, probaću pa ću se javiti.

sve sam probao sem da stavim uslove pre $\theta$ petlje

ovo ostalo sam lupao jer znam da treba da driftuje ka 0, ali nikako nece, pa sam probao svakakve gluposti.

a, šta da radim sa ovim $\lim_{\theta \to 0} {\frac {1}{\theta}\frac{\partial}{\partial\theta}(\theta\frac{\partial p(\theta,z)}{\partial\theta}) = 2\frac{\partial^2p(\theta,z)}{\partial^2\theta} _{za\theta=0} }$

Re: Metod konačnih razlika u Fortranu 24.02.2010. u 14:14 · pre 203 meseca #7

Pardon, treba abs(i * korakuglovni - ugao) .lt. korakuglovni, pošto je Dtheta u radijanima.

Citat:

bdrljaca: a, šta da radim sa ovim [...]

Pa ništa. Jeste singularitet u $\theta = 0$ , ali se automatski izbegava realizacijom Nojmanovog uslova $p_{0,j} = p_{1,j}$ .

Re: Metod konačnih razlika u Fortranu 24.02.2010. u 19:19 · pre 203 meseca #8

Da, i meni je bilo jasno da se duplira uslov, ali pošto prvi put radim, rekoh da proverim. nisam stigao još da proverim ali ću se javiti ako bude problema.
Pozdrav

Re: Metod konačnih razlika u Fortranu 25.02.2010. u 10:26 · pre 203 meseca #9

Isprobao sam i radi. Hvala ti puno, a ako budem negde zakočio javiću se.
Treba da vidim kako da izračunam inverznu Furijeovu u slučaju kad mi je $p(\theta,z,t)$ , tj. iz $p(\theta,z,t)$ prelazim u $p(\theta,z,\omega)$ pomoću Furijeove, pa prethodno pomenutu numeriku koristim za računanje u frekventnom domenu (modifikovanu naravno) i onda treba da se vratim u vremenski domen.
Meni jedino pada da fitujem funkciju koju sam dobio, pa da analitički tražim inverznu. Jel ti pada na pamet nešto elegantnije?

Re: Metod konačnih razlika u Fortranu 26.02.2010. u 11:20 · pre 203 meseca #10

Nisam ništa radio sa Furijeovim transformacijama, ali kao što postoji diskretna Furijeova transformacija, mora da postoji i inverzna diskretna Furijeova transformacija. Ona bi načelno bila bezbednija (u smislu da se ne unose spoljašnji elementi sumnjive prihvatljivosti) od analitičke aproksimacije pa transformacije, a možda i lakša za izvedbu (manje koda).

Re: Metod konačnih razlika u Fortranu 26.02.2010. u 12:09 · pre 203 meseca #11

Znam, nego nikad nisam radio :).
Još nešto da te pitam, da ne znaš kako bi izgledalo u Fortranu kada bih umesto Delta funkcije na ulazu ubacio Gausijan?

Forum kao aplikacija na telefonu

Srodne teme