Finding the maximum value of x in a differential equation.

[theBEAST]

Homework Statement

I want to find the maximum value of x given that the differential equation relating x is:

x' = 0.7548 - 599.49x^6 -2.4038x^2 - 0.12236x^1.5

Where x = 0 at t = 0.

I have to do this with MATLAB using Euler's Method and here is my code:

x = 0
a = 0 % initial t value
b = 1000 % final t value
N = 1000 % number of intervals

for i = 1:N-1
    x(i+1) = x(i) + (b-a)/N * (0.7546-599.49*x(i)^6 - 2.4038*x(i)^2 - 0.12236*x(i)^1.5)
end

max(x)

I get an imaginary number for the value... And most of the numbers in the vector x are NaN + NaNi.

Does anyone know what I am doing wrong?

 

[haruspex]

Looks like a, b and N never change, so the only use you make of them in the code reduces to the constant 1.

 

[theBEAST]

Looks like a, b and N never change, so the only use you make of them in the code reduces to the constant 1.

Yeah so shouldn't the code output the value of x at t = (0,1,2,...,1000)?

 

[haruspex]

Yeah so shouldn't the code output the value of x at t = (0,1,2,...,1000)?

If that's what you want, yes, but how did you pick that range? Your algorithm is shooting way beyond where you want to be and dives down to where x goes negative. Once it has done that the power of 1.5 will produce imaginary results.
Can you think of a way of determining a smaller range that should include a solution?

 

[theBEAST]

If that's what you want, yes, but how did you pick that range? Your algorithm is shooting way beyond where you want to be and dives down to where x goes negative. Once it has done that the power of 1.5 will produce imaginary results.
Can you think of a way of determining a smaller range that should include a solution?

Alright, but before figuring out a smaller range... I am having problems with my code already.

So if I use a = 0 and b = 1000 I get a vector x where x(2) = 0.7546. In other words x value at t = 1 gives x = 0.7546.

NOW if I decide to use a = 0 and b = 1. Then I should expect that max(x) gives me some value greater than or equal to 0.7546. Instead it gives me 0.3074...

From this I think there is something wrong with my code...?

 

[haruspex]

So if I use a = 0 and b = 1000 I get a vector x where x(2) = 0.7546. In other words x value at t = 1 gives x = 0.7546.

NOW if I decide to use a = 0 and b = 1. Then I should expect that max(x) gives me some value greater than or equal to 0.7546. Instead it gives me 0.3074...

With b=1000, you are effectively considering t up to 1000. x[2] corresponds to t=2. With b = 1, you are only considering t up to 1. x[2] now corresponds to t=2/1000.
Even if two setting of b encompass the maximum you're looking for, you would expect different values because of error accumulation, and it's not obvious which will give the greater value.

 

[theBEAST]

With b=1000, you are effectively considering t up to 1000. x[2] corresponds to t=2. With b = 1, you are only considering t up to 1. x[2] now corresponds to t=2/1000.
Even if two setting of b encompass the maximum you're looking for, you would expect different values because of error accumulation, and it's not obvious which will give the greater value.

Ah, so the one where b = 1 would give me a more accurate value but it may not be the interval where the max is located... I tried to figure out the interval by changing the values of b but as you said that would cause error to accumulate. Right now I am not sure what to do...

As I increase b, it oscillates more and more.

 

[haruspex]

Well, I would not be using Euler's method to solve this problem. You want the point where x'=0, so that's the same as solving 0.7548 - 599.49x⁶ -2.4038x² - 0.12236x^1.5 = 0, or, as a polynomial in y = √x, 0.7548 - 599.49y¹² -2.4038y⁴ - 0.12236y³ = 0.
At the least, you could use that to get a handle on the range. Just find some smallish y, y_lim, which makes the polynomial go negative. (It's going to be less than 1, right?) Set b-a = y_lim².

 

[theBEAST]

Well, I would not be using Euler's method to solve this problem. You want the point where x'=0, so that's the same as solving 0.7548 - 599.49x⁶ -2.4038x² - 0.12236x^1.5 = 0, or, as a polynomial in y = √x, 0.7548 - 599.49y¹² -2.4038y⁴ - 0.12236y³ = 0.
At the least, you could use that to get a handle on the range. Just find some smallish y, y_lim, which makes the polynomial go negative. (It's going to be less than 1, right?) Set b-a = y_lim².

LOL wow I can't believe I forgot about that. When the derivative is zero you get either a max or min. In this case it is a max because the second derivative at the solution I got x = 0.30748 is negative!

And this corresponds to the value I get from Euler's!

Thanks!

 

[haruspex]

LOL wow I can't believe I forgot about that. When the derivative is zero you get either a max or min. In this case it is a max because the second derivative at the solution I got x = 0.30748 is negative!

And this corresponds to the value I get from Euler's!

Thanks!

The only reason I didn't mention this before is that you said in the OP that you had to use Euler's method. Was that not true?

 

[theBEAST]

The only reason I didn't mention this before is that you said in the OP that you had to use Euler's method. Was that not true?

Oh, the question said I had to use MATLAB so the first thing that popped up in my mind was either euler's method or ode45. But then after what you said, you can also use MATLAB to solve for the zeros of the equation, which is a much better method.