ODE with diff? DE in x(t), w/ forcing y(t), where y'(t) also appears.

Hi folks, I'm implementing a problem here in Matlab, and I'm having a bit of trouble with the right approach. The bulk of my experience is with Mathematica, where NDSolve is a bit of a black box, but I'm trying this in Matlab as a learning exercise.

I'm looking to solve an ODE over time, with a specific forcing function, for which the derivative also appears. This is for a simple mass-spring-damper with seismic excitation (not force excitation).

m*x''(t) + c*x'(t) + k*x(t) + m*g = k*y(t) + c*y'(t)

Here, y(t) is entirely specified before hand, where y(t) is a generated input that is meant to be general. So, for RHS = 0, I can get this system working just fine with ode45; that's easy enough.

What I'm uncertain of is how to deal with the y(t) and y'(t). Since y(t) is specified ahead of time, should I also calculate a yprime(t) = diff(y)/tstep; , and then inside the odefunc, do an interp1 to look up the value at the given time? This seems a little clunky; if anyone knows of a better way, please let me know.

Answer 1

[Teja Muppirala]

It might seem counterintuitive, but even though y'(t) appears in the ODE, you don't actually need to calculate it. No DIFF, no GRADIENT, no analytical derivatives, no anything. First let me show you that this is true, then I will show you why it is true.

Just as an example, say, m = c = g = k = 1. We'll say y(t) = cos(t), and we'll take some arbitrary initial conditions x(0) = 2 and x'(0) = 3.

% Method 1, your current method, explicitly differentiating to get -sin(t):
F = @(t,x) [0 1; -1 -1]*x + [0;-1 + cos(t) - sin(t)];
x0 = 2; xp0 = 3; IC = [x0; xp0];
[tout, xout1] = ode45(F,[0:0.1:10],IC);
figure, plot(tout,xout1(:,1));
hold all

% Method 2, without differentiating y in any way:
F = @(t,z) [-1 1; -1 0]*z + [1;1]*[cos(t)] + [-t;0];
IC = [x0; x0 + xp0 - cos(0)];
[tout, xout2] = ode45(F,[0:0.1:10],IC);
plot(tout,xout1(:,1),'r:'), legend('Method 1', 'Method 2');

You can see they give exactly the same answer, and it is true in general, you will never have to differentiate y for this sort of linear ODE.

Now, I will will show you how I derived this. Start with your equation:

m*x''(t) + c*x'(t) + k*x(t) + m*g = k*y(t) + c*y'(t)

Step 1. Divide through by m. I will call k/m = K and c/m = C.

x'' + C*x' + K*x + g = K*y + C*y'

Step 2. Rearrange to get x'' by itself.

x'' = -C*x' - K*x + K*y + C*y' - g

Step 3. Integrate this up once with respect to time, assuming zero intial conditions for everything (don't worry, we'll deal with them later). I will call int(y) = Y and int(x) = X.

x' = -C*x - K*X + K*Y + C*y - g*t

Step 4. Define z1 = x and rewrite.

z1' = -C*z1 - K*X + K*Y + C*y - g*t

Step 5. Define -K*X + K*Y = z2, then differentiate z2.

z1' = -C*z1 + z2 + C*y - g*t

z2' = -K*x + K*y = -K*z1 + K*y

Step 6. If you let z = [z1;z2], the final equation is:

z' = [-C 1; -K 0]*z + [C;K]*y + [-g*t; 0]

And, as we defined in step 4, your final answer x(t) is given by the first state, z1.

Step 7. How are z(0) and z'(0) related to x(0) and x'(0)? Figure out how to map initial conditions.

z1(0) is simply equal to x(0).

z2(0) is a bit more tricky. Plug in t=0 into the equation for z1'.

z1'(0) = -C*z1(0) + z2(0) + C*y(0) - g*0 = -C*z1(0) + z2(0) + C*y(0)

Rearrange this to solve for z2(0).

z2(0) = z1'(0) + C*z1(0) -C*y(0)

z2(0) = x'(0) + C*x(0) -C*y(0)

DONE

There you have it. This is a very general method, and it is used in particular to generate state space equations from linear ODEs. See http://en.wikipedia.org/wiki/State_space_representation#Canonical_realizations

Answer 2

[Jan]

There is a small improvement by using gradient(y,t) instead of diff(x)/tstep. Then notice, that a linear interpolation is not differentiable, such that there are discontinuities in the function to be integrated. Matlab's integrators cannot handle this correctly, see http://www.mathworks.com/matlabcentral/answers/59582#answer_72047. So try a cubic interpolation instead.

If you determine y(t), is there a chance to find its derivative also?

Answer 3

[Andrew]

Thanks for the answer, Jan. gradient() makes a lot more sense than the way I was doing it with diff, since I had to pad the array, since diff returns a length-1 array.

I also switched the interpolation over to cubic interpolation; it was actually working OK with just linear interpolation, as far as I know. That might not be the case once I start applying more complicated dynamics.

A detail that I didn't include in full in my original question was: What is the "right" way to implement this forcing function? I ended up passing the array of time points, and all of y(t) and the re-calculated y'(t) as arguments into the odefunc. That seems clumsy, but it was the only way I could think of doing it. If there is another, cleaner way of doing that, I'd love to know.

For reference, here's the code that I currently have working: Problem setup: http://pastebin.com/7gS4N3te ODE function: http://pastebin.com/Nbx58nZ1