SHOW:
|
|
- or go back to the newest paste.
| 1 | module params | |
| 2 | real:: L, K , M, g | |
| 3 | real :: a, b, relerr, abserr, result, errest, flag !Starting values for QUANC8 | |
| 4 | real :: T, TOUT, Y(4), R_WORK(28) !For RKF45 | |
| 5 | integer :: IFLAG, IWORK(5), NEQN !For RKF45 | |
| 6 | real :: U_STATIC(7), U_DYNAMIC(7), fm !Vectors for U and their quadratic sum | |
| 7 | DATA U_STATIC/0.0, 0.303, -0.465, 0.592, -0.409, 0.164, 0.180/ !Input data for U | |
| 8 | real :: fm_a , fm_b , tol !For fmin | |
| 9 | ||
| 10 | end module params | |
| 11 | ||
| 12 | program kurs15 | |
| 13 | use params | |
| 14 | implicit none | |
| 15 | integer i, nofun; | |
| 16 | real:: fmin | |
| 17 | ||
| 18 | interface | |
| 19 | ||
| 20 | real function for_l(x) | |
| 21 | real x | |
| 22 | end function for_l | |
| 23 | ||
| 24 | subroutine diff(T1, Y1, YP) | |
| 25 | real T1, Y1(4), YP(4) | |
| 26 | end subroutine diff | |
| 27 | ||
| 28 | real function sub_fmin(x) | |
| 29 | real x | |
| 30 | end function sub_fmin | |
| 31 | ||
| 32 | end interface | |
| 33 | ||
| 34 | M = 1.0 | |
| 35 | g = 9.81 | |
| 36 | fm_a = 36.0 | |
| 37 | fm_b = 46.0 | |
| 38 | tol = 1.0E-04 | |
| 39 | a = 0.0 | |
| 40 | b = 1.0 | |
| 41 | relerr = 1.E-05 | |
| 42 | abserr = 1.E-05 | |
| 43 | NEQN = 4 | |
| 44 | ||
| 45 | call QUANC8(for_l, a, b, abserr, relerr, result, errest, nofun, flag) | |
| 46 | L = 0.6836896 * result | |
| 47 | write(*,*) "L = ", L | |
| 48 | ||
| 49 | K = FMIN(fm_a, fm_b, sub_fmin, tol) | |
| 50 | write(*,*) "FMIN = ", K | |
| 51 | ||
| 52 | T = 0.0 | |
| 53 | TOUT= 0.0 | |
| 54 | IFLAG = 1 | |
| 55 | ||
| 56 | Y(1) = 0.0 | |
| 57 | Y(2) = 0.0 | |
| 58 | Y(3) = 0.0 | |
| 59 | Y(4) = 4.0 | |
| 60 | fm = 0.0 | |
| 61 | do i=1, 7 | |
| 62 | CALL RKF45(diff,NEQN,Y,T,TOUT,RELERR, ABSERR, IFLAG, R_WORK, IWORK) | |
| 63 | ||
| 64 | write(*,*) "T = ", T, "IFLAG = ", IFLAG | |
| 65 | write(*,*) "Y = ", Y | |
| 66 | write(*,*) "----------------------------------" | |
| 67 | TOUT = T + 0.4 | |
| 68 | fm = fm + (Y(1) - U_STATIC(i))**2 | |
| 69 | end do | |
| 70 | write(*,*) "fm = ", sqrt(fm)/7 | |
| 71 | ||
| 72 | !------------------------------------------------------------------- | |
| 73 | ||
| 74 | end program kurs15 | |
| 75 | ||
| 76 | real function for_l(x) | |
| 77 | implicit none | |
| 78 | real x | |
| 79 | for_l = exp(x) | |
| 80 | end function for_l | |
| 81 | ||
| 82 | subroutine diff(T1, Y1, YP) | |
| 83 | ||
| 84 | use params | |
| 85 | implicit none | |
| 86 | real T1, Y1(4), YP(4) | |
| 87 | real xd, yd, od, zd | |
| 88 | xd = Y1(1) | |
| 89 | od = Y1(2) | |
| 90 | yd = Y1(3) | |
| 91 | zd = Y1(4) | |
| 92 | YP(1) = yd | |
| 93 | YP(2) = zd | |
| 94 | YP(3) = -K/M*xd - g*(1 - cos(od)) + (L+xd)*zd**2 | |
| 95 | YP(4) = -g/(L+xd)*sin(od) - 2/(L+xd)*yd*zd | |
| 96 | end subroutine diff | |
| 97 | ||
| 98 | real function sub_fmin(x) | |
| 99 | use params | |
| 100 | ||
| 101 | real x | |
| 102 | interface | |
| 103 | subroutine diff(T1, Y1, YP) | |
| 104 | real T1, Y1(3), YP(3) | |
| 105 | end subroutine diff | |
| 106 | end interface | |
| 107 | K = x | |
| 108 | fm = 0.0; | |
| 109 | T = 0.0 | |
| 110 | TOUT= 0.0 | |
| 111 | IFLAG = 1 | |
| 112 | ||
| 113 | Y(1) = 0.0 | |
| 114 | Y(2) = 0.0 | |
| 115 | Y(3) = 0.0 | |
| 116 | Y(4) = 4.0 | |
| 117 | do i=1, 7 | |
| 118 | CALL RKF45(diff,NEQN,Y,T,TOUT,RELERR, ABSERR, IFLAG, R_WORK, IWORK) | |
| 119 | TOUT = T + 0.4 | |
| 120 | fm = fm + (Y(1) - U_STATIC(i))**2 | |
| 121 | end do | |
| 122 | sub_fmin = fm | |
| 123 | end function sub_fmin |
