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 |