1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Fortran77 woes- single parameter least squares minimisation stuck!

  1. Jan 26, 2014 #1
    Hey guys! Having issues with a code I'm writing. It doesn't fit the typical format and I may have already accidentally posted it in the wrong forum, so here's a copy of what was there with the newest version of the code. It'll (eventually) be a four-parameter least squares minimisation, but for the moment I'm just trying to get the single parameter version to work. I've got the first module done but, when I run it it gets caught up in a loop and doesn't actually minimise anything! Good start huh?

    Here's my code, bless me with your insight!

    PROGRAM assignment

    ! A program designed to fit experiemental data, using the method
    ! of least squares to minimise the associated chi-squared and
    ! obtain the four control parameters A,B,h1 and h2.
    !*****************************************************************

    IMPLICIT NONE
    INTEGER i,t0
    DOUBLE PRECISION t(17),t2(17),Ct(17),eCt(17),Ctdiff(17),c(17)
    DOUBLE PRECISION Aa,Ab,Ba,Bb,h1a,h1b,h2a,h2b,chisqa,chisqb,dchisq
    DOUBLE PRECISION deltah,Cs(17)

    OPEN(21, FILE='data.txt', FORM='FORMATTED', STATUS='OLD')
    DO i=1,17
    READ(21,*)t(i),Ct(i),eCt(i)
    END DO
    CLOSE(21)

    !Read in data.txt as three one dimensional arrays.
    !*****************************************************************
    !OPEN(21, FILE='outtest.txt', FORM='FORMATTED', STATUS='NEW')
    !DO i=1,17
    ! WRITE(21,*)t(i),Ct(i),eCt(i)
    !END DO
    !CLOSE(21)
    !
    !Just to check input file is being read correctly.
    !*****************************************************************
    !****************************Parameters***************************

    Aa= 0
    Ba= 0
    h1a= 0
    h2a= 0

    !*****************************************************************
    !**********************Minimising Lamda1 (h1)*********************
    deltah= 0.0001
    h1b= 0.001
    h1a= 0
    DO 10
    chisqa= 0
    DO 20 i= 1,17
    Cs(i)= exp(-h1a*t(i))!*Aa !+ Ba*exp(-h2a*t(i))
    Ctdiff= Ct - Cs
    c= Ctdiff**2/eCt**2
    chisqa= chisqa + c(i)
    20 END DO
    !Use initial h1 value of 0 to calculate start-point chisq.

    chisqb= 0
    DO 30 i= 1,17
    h1b= h1b + deltah
    Cs(i)= exp(-h1b*t(i))!*Ab !+ Bb*exp(-h2b*t(i))
    Ctdiff(i)= Ct(i) - Cs(i)
    c= Ctdiff(i)**2/eCt(i)**2
    chisqb= chisqb + c(i)
    30 END DO
    !First-step h1 used to find competing chisq for comparison.
    WRITE(6,*) 'Chi-squared a=', chisqa,'for Lamda1=',h1a
    WRITE(6,*) 'Chi-squared b=', chisqb,'for Lamda1=',h1b
    !Prints the two calculated chisq values to screen.
    dchisq= chisqa - chisqb
    IF (dchisq.GT.0) THEN
    h1a= h1b
    ELSE IF (dchisq.LE.0) THEN
    deltah= deltah-((deltah*2)/100)

    END IF
    IF (chisqb.LE.6618.681) EXIT
    10 END DO
    WRITE(6,*) 'Chi-squared is', chisqb,'for Lamda1=',h1b
    !*****************************************************************

    END PROGRAM assignment


    Thanks for the help.
    Nathan

    EDIT: I should be getting a chi-squared of 6618.681 from this, but it's just stuck between 6921.866 and 6920.031. Help!
     
  2. jcsd
  3. Jan 26, 2014 #2

    Mark44

    Staff: Mentor

    I added some comments where I think you're going wrong. My comments look like this: !!! ***
    In one comment, you have three assignment statements where you are using array variables, but without indexes. This is wrong. Everywhere you have an array variable in a calculation, it should have an index. E.g., it should be Ct(i), not just Ct.

    In my other comments, you have a line of code that might be wrong. It looks like you commented it out, but I'm not sure what you did.

    Your variable names are terrible! For example, Aa,Ab,Ba,Bb,h1a,h1b. Can't you come up with some names that are self-explanatory? Also, why does your code meander all across the page? This makes it harder to see which statements are inside IF blocks and DO loops and such.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Fortran77 woes- single parameter least squares minimisation stuck!
Loading...