SHOW:
|
|
- or go back to the newest paste.
| 1 | PROGRAM awp | |
| 2 | ||
| 3 | IMPLICIT NONE | |
| 4 | ||
| 5 | INTEGER:: N, i, j, idat | |
| 6 | DOUBLE PRECISION:: t0, tend, E_1, E_2, dt_1, dt_2, o,dy(4),y(4),t,yex(2),v0,r0,T_0,pi | |
| 7 | CHARACTER::datei | |
| 8 | pi=acos(-1.) | |
| 9 | ||
| 10 | t0=0. | |
| 11 | tend=0.6 | |
| 12 | v0=10. | |
| 13 | r0=1. | |
| 14 | T_0=2.*pi*r0*1./v0 | |
| 15 | ||
| 16 | ||
| 17 | OPEN(UNIT=100,FILE='ausgabe_ec.csv',STATUS='UNKNOWN') | |
| 18 | WRITE(100,'(A)') "Streifen; analytisch; numerisch; Schrittweite; Fehler; Ordnung" | |
| 19 | ||
| 20 | idat=1 | |
| 21 | DO i=4,10 | |
| 22 | N=2**i | |
| 23 | dt_2=DBLE(tend-t0)/DBLE(N) | |
| 24 | y=(/ r0,0.,0.,v0 /) | |
| 25 | t=0. | |
| 26 | E_2=0. | |
| 27 | ||
| 28 | WRITE(datei,'(I1)')idat | |
| 29 | idat=idat+1 | |
| 30 | write(*,*)idat | |
| 31 | OPEN(UNIT=101,FILE='ec_erg'//TRIM(datei)//'.csv',STATUS='UNKNOWN') | |
| 32 | !OPEN(UNIT=101,FILE='rk_erg'//TRIM(datei)//'.csv',STATUS='UNKNOWN') | |
| 33 | WRITE(101,*)"Zeit;x;xex;y;yex;vx;vy" | |
| 34 | ||
| 35 | !Anfangswerte | |
| 36 | yex=yexakt(t,r0,T_0) | |
| 37 | WRITE(101,'(F12.8,A,F12.8,A,F12.8,A,F12.8,A,F12.8,A,F12.8,A,F12.8)') & | |
| 38 | & t,";",y(1),";",yex(1),";",y(2),";",yex(2),";",y(3),";",y(4) | |
| 39 | ||
| 40 | DO j=0,N-1 | |
| 41 | dy = eulercauchy(y,dt_2,t) | |
| 42 | !dy = rungekutta(y,dt,t) | |
| 43 | y = y + dy | |
| 44 | t = t + dt_2 | |
| 45 | yex=yexakt(t,r0,T_0) | |
| 46 | WRITE(101,'(F12.8,A,F12.8,A,F12.8,A,F12.8,A,F12.8,A,F12.8,A,F12.8)') & | |
| 47 | & t,";",y(1),";",yex(1),";",y(2),";",yex(2),";",y(3),";",y(4) | |
| 48 | ||
| 49 | E_2 = E_2 + (y(1)-yex(1))**2 | |
| 50 | END DO | |
| 51 | CLOSE(101) | |
| 52 | ||
| 53 | E_2=sqrt(E_2/N) | |
| 54 | ||
| 55 | IF (i>4) THEN | |
| 56 | o=DLOG(E_1/E_2)/DLOG(dt_1/dt_2) | |
| 57 | ELSE | |
| 58 | o=0. | |
| 59 | END IF | |
| 60 | ||
| 61 | WRITE(100,'(I8, A, F12.8, A, F12.8, A, F12.8, A, E14.7, A, F10.3)')N,& | |
| 62 | &";", yex(1), ";", y(1), ";", dt_2, ";", E_2, ";", o | |
| 63 | ||
| 64 | E_1=E_2 | |
| 65 | dt_1=dt_2 | |
| 66 | END DO | |
| 67 | ||
| 68 | CLOSE(100) | |
| 69 | ||
| 70 | write(*,*) y(1) | |
| 71 | write(*,*) y(2) | |
| 72 | write(*,*) y(3) | |
| 73 | write(*,*) y(4) | |
| 74 | ||
| 75 | CONTAINS | |
| 76 | ||
| 77 | ||
| 78 | ||
| 79 | ||
| 80 | ||
| 81 | function f(y,t) | |
| 82 | real::f(4) | |
| 83 | real::B | |
| 84 | real, intent(in)::y(4),t | |
| 85 | B=-10. | |
| 86 | write(*,*) "vx,vy: ",y(3)," ",y(4) | |
| 87 | f=(/ y(3),y(4),y(4)*B,-y(3)*B /) | |
| 88 | end function f | |
| 89 | ||
| 90 | function yexakt(t,r,T_0) | |
| 91 | real::yexakt(2),pi | |
| 92 | real, intent(in)::t,r,T_0 | |
| 93 | pi=acos(-1.) | |
| 94 | yexakt=(/ r*cos(2.*pi/T_0*t), r*sin(2.*pi/T_0*t) /) | |
| 95 | end function yexakt | |
| 96 | ||
| 97 | ||
| 98 | function eulercauchy(y,dt,t) | |
| 99 | real::eulercauchy(4) | |
| 100 | real, intent(in) :: y(4), dt, t | |
| 101 | eulercauchy=dt*f(y,t) | |
| 102 | end function eulercauchy | |
| 103 | ||
| 104 | function rungekutta(y,dt,t) | |
| 105 | real::rungekutta(4) | |
| 106 | real, intent(in) :: y(4), dt, t | |
| 107 | real:: k1(4), k2(4), k3(4), k4(4) | |
| 108 | ||
| 109 | k1=f(y,t) | |
| 110 | k2=f(y+dt/2*k1,t+dt/2.) | |
| 111 | k3=f(y+dt/2*k2,t+dt/2.) | |
| 112 | k4=f(y+dt*k3,t+dt) | |
| 113 | rungekutta=dt/6.*(k1+2.*(k2+k3)+k4) | |
| 114 | ||
| 115 | end function rungekutta | |
| 116 | ||
| 117 | ||
| 118 | END PROGRAM |
