MATLAB: How to force a non-trivial solution for a second order ODE

MATLABnon-trivialrunge kuttasecond order ode

I am trying to solve the second order ODE x''=[ 9/4-(12/(5-3exp(2t)))exp(2t)] x(t), with ICs x'=x=0 at t=-1 for -1<=t<=0, using Runge-Kutta 4. Here is my MWE:

a = -1; b = 0; N = 1000;

t = zeros(1, N);

w1 = zeros(1, N);

w2 = zeros(1, N);

w1(1) = 0; w2(1) = 0;

h = (b – a)/N;

t(1) = a;

F1 = @(t, x1, x2) x2;

F2 = @(t, x1, x2) ((12*exp(2*t))/(3*exp(2*t) – 5) + 9/4)*x1;

for i = 1:(N-1)

K11 = h*(F1(t(i), w1(i), w2(i)));

K21 = h*(F2(t(i), w1(i), w2(i)));

K12 = h*(F1(t(i) + 0.5*h, w1(i) + 0.5*K11, w2(i)+ 0.5*K21));

K22 = h*(F2(t(i) + 0.5*h, w1(i) + 0.5*K11, w2(i)+ 0.5*K21));

K13 = h*(F1(t(i) + 0.5*h, w1(i) + 0.5*K12, w2(i)+ 0.5*K22));

K23 = h*(F2(t(i) + 0.5*h, w1(i) + 0.5*K12, w2(i)+ 0.5*K22));

K14 = h*(F1(t(i) + h, w1(i) + K13, w2(i)+ K23));

K24 = h*(F2(t(i) + h, w1(i) + K13, w2(i)+ K23));

w1(i+1) = w1(i) + (K11 + 2*K12 + 2*K13 + K14)/6;

w2(i+1) = w2(i) + (K21 + 2*K22 + 2*K23 + K24)/6;

t(i+1) = a + i*h;

end

My issue is that the trivial solution x=0 is returned. I have tried the code with various other second order ODEs that I know the solution for and I get the expected results. If I make the ICs tiny instead of zero (say 0.000001) it seems to force a correct-ish looking, but obviously incorrect, solution. Any help is appreciated. I am new to Matlab. I get the same thing using ode45. Perhaps there is no non-trivial solution but there should be.

Best Answer

Your equation is x''=[ 9/4-(12/(5-3exp(2t)))exp(2t)] x(t) with initial conditions x'=x=0. You are a racecar at the starting line (since x = 0), the car is not moving (since x' = 0 its speed is 0 miles [or kilometers] per hour), and your foot is off the accelerator (evaluating the right hand side gives x'' = 0 regardless of what value t has.)

Is your car ever going to move under those circumstances? It will when they bring the tow truck over to pull you out of the way so others can race, but that's outside the scope of this ODE.

As for your simpler example, with x'' = -x+6*t+t^3, that's different than this example. You're not multiplying by x.

A closer example would be x'' = (-x+6*t+t^3)*x and again, with the initial conditions x = 0 and x' = 0 you never start accelerating.

Related Solutions

MATLAB: Runge kutta 4 with two ODE’s – function inside a function

Since W is a function of time, it needs to match the time that you are using in each particular line of code. E.g., take this line:

K11 = fy1(t(i)      ,y1(i)        , y2(i), W(i));

The time for this line is t(i), so the W will be fishing_load_factor(t(i)). E.g.,

K11 = fy1(t(i)      ,y1(i)        , y2(i), fishing_load_factor(t(i)));

For the above lines you get the same result, but now look at this line:

K21 = fy1(t(i)+h/2  ,y1(i)+h/2*K11, y2(i)+h/2*K11, W(i));

Here the time is different, but you didn't adjust your W accordingly. It should be this instead:

K21 = fy1(t(i)+h/2  ,y1(i)+h/2*K11, y2(i)+h/2*K11, fishing_load_factor(t(i)+h/2));

Similarly for the other lines of code involving W.

In addition, the times of these lines are incorrect:

    K41 = fy1(t(i)+2    ,y1(i)+h*K31   , y2(i)+h*K31, W(i));
    K42 = fy1(t(i)+2    ,y1(i)+h*K32  , y2(i)+h*K32);

The times should be t(i)+h instead of t(i)+2.

Then, the factors for the K31 and K32 terms in these lines are incorrect:

    y1(1+i) = y1(i)+h/6*(K11 + 2*K21 + K21*K31 + K41);
    y2(1+i) = y2(i)+h/6*(K12 + 2*K22 + K22*K32 + K42);

Instead of K21*K31 and K22*K32 they should be 2*K31 and 2*K32 respectively.

FInally, you are using the derivative calculations with the wrong states. y1 derivatives should be calculated with the fy1 function and applied to the y1 state, and y2 derivatives should be calculated with the fy2 function and applied to the y2 state. But look what you have done with these lines:

    K21 = fy1(t(i)+h/2  ,y1(i)+h/2*K11, y2(i)+h/2*K11, W(i));
    K22 = fy1(t(i)+h/2  ,y1(i)+h/2*K12, y2(i)+h/2*K12);

K11 came from a fy1 calculation, so it should only be used with the y1 state, but you are using it with the y2 state in this expression:

    y2(i)+h/2*K11

That expression should be using the derivative that came from the fy2 function, namely K12:

    y2(i)+h/2*K12

And then your K22 calculation, which will be used with the y2 state, uses the fy1 function instead of the fy2 function. These two lines should be this instead:

    K21 = fy1(t(i)+h/2  ,y1(i)+h/2*K11, y2(i)+h/2*K12, fishing_load_factor(t(i)+h/2));
    K22 = fy2(t(i)+h/2  ,y1(i)+h/2*K11, y2(i)+h/2*K12);

So, LOTS of these types of errors in your code that need to get corrected. Just follow this rule:

y1 derivatives always use the fy1 function and y1 always gets propagated with the Kx1 derivatives

y2 derivatives always use the fy2 function and y2 always gets propagated with the Kx2 derivatives

I would add that it looks like your h is way too big to pick up the change points in your W(t) function properly. You will need to shrink h by a fair amount to make sure you pick up those W(t) change points adequately.

MATLAB: Error using plot Vectors must be the same length.

Some comments

1) Q is defined twice as symbolic variable, but it is actually an array

2) in the first loop N = 4, so length(t) = 4 and length(E1) = 4; [btw: the real length is 5 because i = 1 is not addressed: your index is always i+1]. In the second loop N = 40, so length(t) = 40 and length(E2) = 40. In the end you use the same t for plottin either E1 and E2.

My suggestion is to use different vectors for the first and second loop, maybe t1 and t2. Then preallocate always your vectors. Finally in the loop, address Q(i) and E1(i) (and E2) instead of Q(i+1) and E1(i+1)

Best Answer

Related Solutions

MATLAB: Runge kutta 4 with two ODE’s – function inside a function

MATLAB: Error using plot Vectors must be the same length.

Related Question