ode45 script not changing initial values (for y) in a plot

I been attempting to graph a system of diff. equations for some time now. I finally got the plots to work, but the initial values are not changing so the graphs are flat.

%create function scripts in separate .m files:
function D = Dtemp(T)

D = (.00000973/298.15^1.75)*T^1.75

function n = Ntemp(T)

n = Ptemp(T)/(0.0138*T)

function p = Ptemp(T)
p=10^(3.93068-(1337.716/(T-82.648)))*(10^5)

%Create the .m file for the Diff. EQ:
function dy = odefun(t,y)
dy = zeros(2,1)
dy(1) = -Dtemp(y(2))*1.21e-25*(Ntemp(y(2))/y(1)*948)
dy(2) = 3*5.78e5*((-Dtemp(y(2))*1.21e-25*(Ntemp(y(2))/y(1)*948))/2000.68*y(1)) + (3*(0.0454/2000.68*y(1))*((573.15-y(2))/948*y(1)))

%Open a new script and run the DE, attempting to plot:

[t,y] = ode45(@odefun,[0 .3869],[1e-6; 298.15])

plot(t,y(:,1),'-o')
figure
plot(t,y(:,2),'-o')

Results:

t =

           0
      0.0097
      0.0193
      0.0290
      0.0387
      0.0484
      0.0580
      0.0677
      0.0774
      0.0871
      0.0967
      0.1064
      0.1161

    y =

        0.0000  298.1500
        0.0000  298.1500
        0.0000  298.1500
        0.0000  298.1500
        0.0000  298.1500
        0.0000  298.1500
        0.0000  298.1500
        0.0000  298.1500
        0.0000  298.1500
        0.0000  298.1500
        0.0000  298.1500
        0.0000  298.1500
        0.0000  298.1500
        0.0000  298.1500
        0.0000  298.1500
        0.0000  298.1500
        0.0000  298.1500
        0.0000  298.1500
        0.0000  298.1500

As you can see my t values are changing, but the y-arrays are not. I'm not sure what the issue is Here. Any help would be much appreciated!

Answer 1

[Star Strider]

I suspect the problem is traceable to the factor ‘1.21e-25’ present in both the ‘dy(1)’ and ‘dy(2)’ calculations. This is essentially zero, and I don’t see anything in your functions that would scale to the same magnitude to cancel them (for example in a multiplication or division with a constant of a similar or inverse magnitude). These are derivatives, and derivatives that small would not change within the range of machine precision (given by the eps function) between iterations.

Since I’m no certain what you’re doing (and it may not be within my areas of expertise anyway), I can’t suggest a work-around.

[Ryan Compton]

Thank you for the reply,

I added the format long but it still looks as though the values really are not iterating.

0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000

If it is in fact the case that the steps are very small, how could I get this work? I also attempted dividing the time steps into 800000 but the graphs were still the same flatlines.

[Star Strider]

The problem is this (in the final odefun call):

dy =
     -1.43427166128484e-19
      1.97479081953014e-17

and since ‘y(1)’ begins at 1E-6, it changes (although slowly), while for ‘y(2)’ that begins with 298.15, the change is so far below the minimum precision in floating-point calculations with respect to the initial value that it actually never does change:

y =
                     1e-06                    298.15
      9.99999999999998e-07                    298.15
      9.99999999999997e-07                    298.15
      9.99999999999996e-07                    298.15

You would have to use some sort of extended-precision arithmetic to see ‘y(2)’ change at all.

John D’Errico submitted an extended-precision arithmetic package to the File Exchange, that requires the Symbolic Math Toolbox, Mapping Toolbox, and Phased Array System Toolbox. I don’t have the last two, so I’ve no experience with it.

Answer 2

[Steven Lord]

The elements of the y array are changing. They're just changing by a very, very small amount. Change the display formatting:

format longg

and evaluate your odefun for t = 0 and y = your initial conditions.

dy = odefun(0, [1e-6; 298.15]);

The elements of dy are not 0 but they are very, very small. The factors of 1.21e-25 in the lines of code that compute the elements of dy and the constant roughly equal to 4.5e-10 in your Dtemp function contribute to making the elements of dy so small.

[Ryan Compton]

Thank you for the reply,

I added the format long but it stil looks as though the values really are not iterating.

0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000
   0.000000010000000   2.981500000000000

If it is in fact the case that the steps are very small, how could I get this work? I also attempted dividing the time steps into 800000 but the graphs were still the same flatlines.