solving odes that matrices should be update in every step.

Hi.I have a question that may be so elementary but so far I haven't found a solution. please if it is possible guide me.

i want to solve the 2nd order diffrential equation for vibration of 3DOF system.the equation is : [m]*{x"} + [C]*{x'} + [K]*{X} = -[M]*{r} d2g/dt2

that M is mass matrix . K , C are stiffness and damping matrices that should be update in every time step.also {r} = [ 0;0;1] and g" is the time-dependent external acceleration.(earthquake records).

for updating K and C matrices,i wrote a cod in matlab contains a for-loop that in every time step (=0.5 second) first update [K] and [C] matrices and then use ode45 to solve equation for time span=0.5 second.but i get suspicious result.for example ode45 output for every time span=0.5 second outputs a 61*1 vector but i think output should be a 6*1 vector like this : [x1; x'1; x2; x'2 ;x3;x'3]

can you tell me where is the problem? maybe i use ode45 in wrong method.

i attach my code to this question.

thank you for your help.

Answer 1

[darova]

I made some changes in your code. Pay attention to refreshing parameters:

function sdot = ODEFUN(t,x)        
    g1 = interp1(TIME,G,t);
    ag1 = interp1(TIME,AG,t);
    %% ...
end

Also i didn't refresh C and K matrices (don't know how)

See attached file

[masoud azar]

Thank you very much darova.your changes were very useful for me to understand how ode45 is working.

but i still dont find a way to refresh matrices C and K.

[darova]

What do they depend on? Do you have a condition/parameter?

[masoud azar]

thank you again darova.

K and C matrices contain Kb and Cb that are stiffness and damping values of base isolator.Kb and Cb depend on base isolator displacement .( Cb=f(Xb) , Kb=g(Xb) ) . at end of any time step Xb is changed and so Cb and Kb should be update for next step.the question is when i using ode45 how to update K and C matrices in every step?

it is necessary to say that temporary i supposed Cb and Kb as constant values. (see line 43 and 44 in your mfile) and after i find the answer of my question they will be corrected as a functions of displacements (Xb).

[Walter Roberson]

.the question is when i using ode45 how to update K and C matrices in every step?

Do not do that. Instead pass in a tspan that corresponds only to the next time interval, passing in "current" Cb and Kb values for that time interval. Take the output of ode45 as the new initial conditions and loop back for the next time interval with updated Cb and Kb values.

Any time you must "update parameters", that implies discontinuities, and ode45 cannot work with discontinuities within any one ode45() call. You need to stop the ode45 at the discontinuity and restart it again.

[masoud azar]

thank you Walter Roberson.

for example suppose i have a earthquack record that contains 10 values of time and acceleration ( t , U"g ) and time step is equal to 0.5 secound.

you means that for 10 values of time and acceleration ( t , U"g ), i should run ode45 for 9 times in loop such that time span for ode45 (tspan) is equal to my time steps(0.5 secound)?

at end of every loop Kb and Cb are updated for new values of Xb .and also output results (displacements and velocities) will be saved as initial conditions (yo) for next loop?

is that right?

[darova]

Just place them in ode function

    function sdot = ODEFUN(t,x)
        Cb = f(x(1));
        Kb = g(x(1))
        C=[Cc,0,0;0,Ci,0;0,0,Cb];
        K=[Kc,0,0;0,Ki,0;0,0,Kb];        
        g1 = interp1(TIME,G,t);
%         ag1 = interp1(TIME,AG,t);
%         F=[mc;mi;mc+mi+mb]*-ag1;
        sdot = zeros(6,1);
        sdot(1:3) = x(4:6);                 % dx/dt (velocity)
        sdot(4:6) = -[0;0;1]*g1 - M1*C*x(4:6) - M1*K*x(1:3);    % dx2/dt2 (acc)
    end

[Walter Roberson]

Just place them in ode function

No!! interp1 defaults to linear interpolation between the points on either side of the probe point. The result is something that is piecewise continuous but whose first derivative holds constant between breakpoints but is discontinuous at each breakpoint. That is not acceptable for ode45 purposes. ode45 uses the tension between quartic (degree 4) interpolation and quintic (degree 5) interpolation to adjust step sizes. quartic and quintic interpolation are not valid for any segment in which the function has discontinuities in the first derivative (or second derivative.)

Therefore, at any place in which a discontinuity is to be expected or can be detected to have occurred, you must stop the ode45() and restart it. Failure to do so is violating the mathematical conditions under which ode45() can be expected to produce a valid answer.

When you using interp1() like that, or you using if inside the calculations, ode45 might or might not notice a problem and complain, but if it does not complain then that does not mean that it produced the right answer.

[Walter Roberson]

I recommend having a good look at the ballode example. It shows using loops that eventually cover the entire time span, using event functions to detect when a discontinuinity has occurred and stopping the run and restarting it.

[masoud azar]

thank you darova and Walter Roberson.

at the moment i am trying to edit and rewrite the matlab code in both methods based on your comments and compare the output results in same inputs.

finally i will report the results here.

i am beginer in matlab and it may take a time.

[darova]

I will be here, waiting