SHARE
TWEET

Untitled

a guest Oct 21st, 2019 83 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. program pre4
  2.  
  3.     implicit none
  4.  
  5.     external p,pol,dpol
  6.  
  7. cºººles matrius estan definides grans per a assegurar que les dades hi caben
  8.  
  9.     double precision p,a,b,h,t,pmat(1:200),x(1:200),dapr(1:200),pol,
  10.  
  11.      +eps,xarrel,a1,b1,xarrel1,dpol,xini(1:5),x0,xarrel3
  12.  
  13.     integer intervals,j,niter,i
  14.  
  15.     common/dades/t
  16.  
  17.  
  18.     open(1,file='P3-19pres.dat')
  19. C apartat 1
  20.  
  21.     write(1,100) '#v',"p'(v)",'p(v)'
  22.     a=1.d0/3.d0+0.1d0
  23.     b=4.d0
  24.     t=0.92d0
  25.     intervals=80   
  26.     h=(a-b)/dble(intervals)
  27.  
  28.     do j=1,intervals
  29.  
  30.       x(j)=a+j*h
  31.  
  32.       pmat(j)=p(x(j))
  33.  
  34.     enddo
  35.  
  36.     call defunt1(intervals,pmat,x,dapr)
  37.  
  38.  
  39.  
  40. cºººaquest bucle nomes escriu a l'arxiu
  41.  
  42.     do j=1,intervals
  43.  
  44.       write(1,101) x(j),dapr(j),pmat(j)
  45.  
  46.     enddo
  47.  
  48.  
  49.  
  50. 100 format(a,15x,a,12x,a)
  51.  
  52. 101 format(e15.8,2x,e15.8,2x,e15.8)
  53.  
  54.  
  55.  
  56. c-------------------------------------------------------------------------
  57.  
  58.  
  59.  
  60. c-------
  61.  
  62.     write(1,*)
  63.  
  64.     write(1,*)
  65.  
  66.  
  67.  
  68.    
  69.  
  70.  
  71.  
  72.  
  73.  
  74.  
  75.  
  76.  
  77.  
  78. cºººsubrutina de derivacio analitica
  79.  
  80.     subroutine defunt1(ndat,funci,x,dfunci)
  81.  
  82.     implicit none
  83.  
  84.     double precision x(ndat),funci(ndat),dfunci(ndat),h
  85.  
  86.     integer ndat,k
  87.  
  88.  
  89.  
  90.     h=x(2)-x(1)
  91.  
  92. cºººels dos extrems, 1 i ndat, corresponents a els limits de l'interval, tenen una formula diferent, els escric fora del bucle
  93.  
  94.     dfunci(1)=(funci(2)-funci(1))/h
  95.  
  96.     dfunci(ndat)=(funci(ndat)-funci(ndat-1))/h
  97.  
  98.     do k=2,ndat-1
  99.  
  100.       dfunci(k)=(funci(k+1)-funci(k-1))/(2.d0*h)
  101.  
  102.     enddo
  103.  
  104.  
  105.  
  106.     end subroutine
  107.  
  108.  
  109.  
  110. c-------------------------------------------------------------------------
  111.  
  112. cºººequacio de van der waals
  113.  
  114.     double precision function p(x)
  115.  
  116.     implicit none
  117.  
  118.     double precision x,t
  119.  
  120.     common/dades/t
  121.  
  122.  
  123.  
  124.     p=8.d0*t/(3.d0*x-1.d0)-3.d0/x**2.d0
  125.  
  126.  
  127.  
  128.     return
  129.  
  130.     end function
  131.  
  132.  
  133.  
  134. c-------------------------------------------------------------------------
  135.  
  136. cºººpolinomi apartat 2
  137.  
  138.     double precision function pol(x)
  139.  
  140.     implicit none
  141.  
  142.     double precision x,t
  143.  
  144.     common/dades/t
  145.  
  146.  
  147.  
  148.     pol=4.d0*t*x**3.d0-9.d0*x**2.d0+6.d0*x-1.d0
  149.  
  150.  
  151.  
  152.     return
  153.  
  154.     end function
  155.  
  156.  
  157.  
  158. cºººderivada del polinomi
  159.  
  160.     double precision function dpol(x)
  161.  
  162.     implicit none
  163.  
  164.     double precision x,t
  165.  
  166.     common/dades/t
  167.  
  168.  
  169.  
  170.     dpol=12.d0*t*x**2.d0-18.d0*x+6.d0
  171.  
  172.  
  173.  
  174.     return
  175.  
  176.     end function
  177.  
  178.  
  179.  
  180. c-------------------------------------------------------------------------
  181.  
  182. cºººbiseccio
  183.  
  184.     subroutine bit1(a,b,eps,niter,fun,xarrel)
  185.  
  186.     implicit none
  187.  
  188.     external fun,dfun
  189.  
  190.     double precision a,b,eps,xarrel,error,xt,fun
  191.  
  192.     integer niter
  193.  
  194.  
  195.  
  196.     error=1000.d0
  197.  
  198.     niter=0
  199.  
  200.  
  201.  
  202. cºººsi no hi ha cap zero en l'interval, termina el proces
  203.  
  204.     if(fun(a)*fun(b).gt.0.d0)then
  205.  
  206.       write(*,*) 'en aquest interval no hi ha cap zero (o hi ha un num
  207.  
  208.      +ero parell de zeros)'
  209.  
  210.       goto 666
  211.  
  212.     endif
  213.  
  214.  
  215.  
  216.     do while(error.gt.eps)
  217.  
  218.       xt=(a+b)/2.d0
  219.  
  220.       niter=niter+1
  221.  
  222. cºººsi f(a) i f(xt) son del mateix signe, posem l'a on abans hi era el xt, 'estretant' el rang
  223.  
  224. cºººsi no son del mateix signe, significa que hi ha un 0 entre ells, i el que era xt ara es b, estretant el rang
  225.  
  226.       if(fun(a)*fun(xt).gt.0.d0)a=xt
  227.  
  228.       if(fun(a)*fun(xt).lt.0.d0)b=xt
  229.  
  230. cºººen el cas de que casualment xt resulta ser un 0, iguala 'a' i 'b' a aquest valor per a que l'error doni 0 i es pari el bucle
  231.  
  232.       if(fun(xt).eq.0) then
  233.  
  234.         a=xt
  235.  
  236.         b=xt
  237.  
  238.       endif
  239.  
  240.       error=dabs(a-b)
  241.  
  242.     enddo
  243.  
  244.     xarrel=a
  245.  
  246.  
  247.  
  248. 666 end subroutine
  249.  
  250.  
  251.  
  252. c-------------------------------------------------------------------------
  253.  
  254. cºººnewton rhapson
  255.  
  256.  
  257.  
  258.     subroutine nrt1(x0,eps,niter,fun,dfun,xarrel)
  259.  
  260.     implicit none
  261.  
  262.     external fun,dfun
  263.  
  264.     double precision x0,eps,xarrel,error,xt,fun,dfun
  265.  
  266.     integer niter
  267.  
  268.  
  269.  
  270.     error=1000.d0
  271.  
  272.     niter=0
  273.  
  274.     do while(error.gt.eps)
  275.  
  276.       xt=x0-fun(x0)/dfun(x0)
  277.  
  278.       niter=niter+1
  279.  
  280.       error=dabs(xt-x0)
  281.  
  282.       x0=xt
  283.  
  284.     enddo
  285.  
  286.     xarrel=x0  
  287.  
  288.  
  289.  
  290.     end subroutine
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top