c----------------------------------------------------------------------- c these are the subroutines for the c c Beam (ODE case) c ODE of dimension 80 c c----------------------------------------------------------------------- subroutine prob(problm,neqn,tbegin,tend,ijac,mljac,mujac) character*(*) problm integer neqn,ijac,mljac,mujac double precision tbegin,tend COMMON /NNNN/ NCOM,NNCOM,NSQ,NQUATR,DELTAS problm = 'beam' neqn = 80 tbegin = 0d0 tend = 5.0d0 ijac = 0 mljac = neqn mujac = neqn C --- SET DEFAULT VALUES N =40 NN=2*N NCOM=N NSQ=N*N NQUATR=NSQ*NSQ NNCOM=NN AN=N DELTAS=1.0D+0/AN return end c----------------------------------------------------------------------- subroutine init(neqn,y) integer neqn double precision y(neqn) integer i do 10 i=1,neqn y(i) = 0d0 10 continue return end c----------------------------------------------------------------------- SUBROUTINE feval(NN,T,TH,DF) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION DF(NN),TH(150),U(150),V(150),W(150) DIMENSION ALPHA(150),BETA(150),STH(150),CTH(150) COMMON /NNNN/ N,NNCOM,NSQ,NQUATR,DELTAS C ----- CALCUL DES TH(I) ET DES SIN ET COS ------------- DO 22 I=2,N THDIFF=TH(I)-TH(I-1) STH(I)=DSIN(THDIFF) 22 CTH(I)=DCOS(THDIFF) C -------- CALCUL DU COTE DROIT DU SYSTEME LINEAIRE ----- IF(T.GT.3.14159265358979324D0)THEN C --------- CASE T GREATER PI ------------ C ---------- I=1 ------------ TERM1=(-3.D0*TH(1)+TH(2))*NQUATR V(1)=TERM1 C -------- I=2,..,N-1 ----------- DO 32 I=2,N-1 TERM1=(TH(I-1)-2.D0*TH(I)+TH(I+1))*NQUATR 32 V(I)=TERM1 C ----------- I=N ------------- TERM1=(TH(N-1)-TH(N))*NQUATR V(N)=TERM1 ELSE C --------- CASE T LESS EQUAL PI ------------ FABS=1.5D0*DSIN(T)*DSIN(T) FX=-FABS FY= FABS C ---------- I=1 ------------ TERM1=(-3.D0*TH(1)+TH(2))*NQUATR TERM2=NSQ*(FY*DCOS(TH(1))-FX*DSIN(TH(1))) V(1)=TERM1+TERM2 C -------- I=2,..,N-1 ----------- DO 34 I=2,N-1 TERM1=(TH(I-1)-2.D0*TH(I)+TH(I+1))*NQUATR TERM2=NSQ*(FY*DCOS(TH(I))-FX*DSIN(TH(I))) 34 V(I)=TERM1+TERM2 C ----------- I=N ------------- TERM1=(TH(N-1)-TH(N))*NQUATR TERM2=NSQ*(FY*DCOS(TH(N))-FX*DSIN(TH(N))) V(N)=TERM1+TERM2 END IF C -------- COMPUTE PRODUCT DV=W ------------ W(1)=STH(2)*V(2) DO 43 I=2,N-1 43 W(I)=-STH(I)*V(I-1)+STH(I+1)*V(I+1) W(N)=-STH(N)*V(N-1) C -------- TERME 3 ----------------- DO 435 I=1,N 435 W(I)=W(I)+TH(N+I)*TH(N+I) C ------- SOLVE SYSTEM CW=W --------- ALPHA(1)=1.D0 DO 44 I=2,N ALPHA(I)=2.D0 44 BETA(I-1)=-CTH(I) ALPHA(N)=3.D0 DO 45 I=N-1,1,-1 Q=BETA(I)/ALPHA(I+1) W(I)=W(I)-W(I+1)*Q 45 ALPHA(I)=ALPHA(I)-BETA(I)*Q W(1)=W(1)/ALPHA(1) DO 46 I=2,N 46 W(I)=(W(I)-BETA(I-1)*W(I-1))/ALPHA(I) C -------- COMPUTE U=CV+DW --------- U(1)=V(1)-CTH(2)*V(2)+STH(2)*W(2) DO 47 I=2,N-1 47 U(I)=2.D0*V(I)-CTH(I)*V(I-1)-CTH(I+1)*V(I+1) & -STH(I)*W(I-1)+STH(I+1)*W(I+1) U(N)=3.D0*V(N)-CTH(N)*V(N-1)-STH(N)*W(N-1) C -------- PUT DERIVATIVES IN RIGHT PLACE ------------- DO 54 I=1,N DF(I)=TH(N+I) 54 DF(N+I)=U(I) RETURN END c----------------------------------------------------------------------- subroutine jeval(neqn,t,y,jac,ldim) integer neqn,ldim double precision t,y(neqn),jac(ldim,neqn) return end c----------------------------------------------------------------------- subroutine solut(neqn,true) integer neqn double precision true(neqn) c c true(1)=-0.005792366591294675 true(2)=-0.016952985507199259 true(3)=-0.027691033129713322 true(4)=-0.038008156558781729 true(5)=-0.047906168597422688 true(6)=-0.057387104352737008 true(7)=-0.066453273134522699 true(8)=-0.075107305819780661 true(9)=-0.083352197654124544 true(10)=-0.091191346546446469 true(11)=-0.098628587001297248 true(12)=-0.105668220037774708 true(13)=-0.112315039540924422 true(14)=-0.1 18574355272698475 true(15)=-0.124452012875526880 true(16)=-0.129954411326390999 true(17)=-0.135088518061004200 true(18)=-0.139861881919410397 true(19)=-0.144282644101482929 true(20)=-0.148359547246256976 true(21)=-0.152101942900106414 true(22)=-0.155519797806080921 true(23)=-0.158623699341992299 true(24)=-0.161424860370167541 true(25)=-0.163935123819275499 true(26)=-0.166166967344037066 true(27)=-0.168133508177817718 true(28)=-0.169848508060189926 true(29)=-0.171326378244038509 true(30)=-0.172582184746215274 true(31)=-0.173631653797526901 true(32)=-0.174491177383960691 true(33)=-0.175177818786287100 true(34)=-0.175709317871242317 true(35)=-0.176104096022807288 true(36)=-0.176381260717507812 true(37)=-0.176560609756417469 true(38)=-0.176662635226010517 true(39)=-0.176708527080694206 true(40)=-0.176720176107510191 true(41)= 0.037473626808570053 true(42)=0.109911788012810762 true(43)=0.179836047447039129 true(44)=0.247242730557127186 true(45)=0.312129382035491301 true(46)=0.374494737701689822 true(47)=0.434338607372647125 true(48)=0.491662035432760524 true(49)=0.546467785483476383 true(50)=0.598760970245279030 true(51)=0.648549361126755851 true(52)=0.695843516905088648 true(53)=0.740657266848912124 true(54)=0.783008174791347177 true(55)=0.822917665884869456 true(56)=0.860411030561688098 true(57)=0.895517550233742218 true(58)=0.928270826293034365 true(59)=0.958708933474210358 true(60)=0.986874782150222219 true(61)=1.012816579967983789 true(62)=1.036587736684594479 true(63)=1.058246826485315033 true(64)=1.077857811432700289 true(65)=1.095490221995530989 true(66)=1.111219164319120026 true(67)=1.125125269269998022 true(68)=1.137294526582397119 true(69)=1.147818025203744592 true(70)=1.156792131966898566 true(71)=1.164318845152484938 true(72)=1.170505992580311363 true(73)=1.175467424328008220 true(74)=1.179323003206967714 true(75)=1.182198586301326345 true(76)=1.184226111211404704 true(77)=1.185543909813440450 true(78)=1.186297084230907673 true(79)=1.186637618874913665 true(80)=1.186724615129383839 return end