solving second order ode problem with ode 45

following ups on this thread https://uk.mathworks.com/matlabcentral/answers/512637-solving-a-mass-spring-damper-system-with-ode45

I am analysing a mass spring damper system too, but mine has multiple degrees of freedom. I have a function that creates a column vector of dydt once time and a state vector are provided and it works fine. if I use it as funcion in ode45 it has problems. Which is the correct form of input vectorial function for ode45?

this is a part of the code I am using with the attacched workspace that has to be loaded by the code to run:

% this is one of many funcions in my code 
% it can solve analysis made by many load steps.
% each load step can be either dynamic or static over a fictitius time
% interval. For now I have a linearised and non-linear g-alpha solvers validated against abaqus, 
% but I wanted to use also a solver with a doubled dof system and use ode45
% though I always struggle to get the right formatting to use embedded
% functions of matlab... and I feel always ashamed, so please help!!
% This is an extrapolated demo of my code and for now it won't
% work for multple dynamic load steps solved with ode45. 
% This demo is thought to help the kind person that is willing to help me and therefore it will stop at the end of te first load steps that should be required
% the problematic part ishighlighted below...

clc
clear all
close all
load('dynamic_explicit_solver_workspace.mat')

t0=0;
%     load_step_counter=1;
    Step{1}=RefConfig;
    d0 = zeros(length(UnknownDispIndex),1);  % initial displacement is zero
    v0 = d0*0;
    a0 = d0*0;
    dMtot=[];
    vMtot=[];
    aMtot=[];
    step_counter=1;% used in the function to moove among different TIME steps.. independently from what type of step if linear or dynamic
    timeline=0; 
    for i=1:input.load_step_number-1 % there are 2 steps in the analysis but I am just interested in getting the first one running 
        % my problem is with even one step of ode45 starting from line 50
        % static load step
        if input.StepSolverType{i}==1 % linear static solver
                %other solver not in this demo
        % dynamic load steps    
        elseif input.StepSolverType{i}==2 % for now only ninear but maybe one day also non linear but on th full model 
                % g_alpha solver
                if input.AlgorithmType{1}==1 
                    if input.reduced{i}==1 
                        err('Galpha can only be used for full models for now')
                    else
                        if input.discreteVScontinuous==1
                            %other solver not in this demo
                            
                        else % discreteVScontinuous==2 ... for continuously varing force vector with time
                            if i==1 % just the very first time step require this update of glob force vec the others are fine as they are updated in the loop
                                Step{step_counter}.GlobForceVec=Step{step_counter}.GlobForceVec{1}(t0+input.deltaT_loadstep{i});
                                Step{step_counter}.GlobForceVecReduced=Step{step_counter}.GlobForceVecReduced{1}(t0+input.deltaT_loadstep{i});
                            else
                            end
                            %other solver not in this demo
                        end
                    end
                    % ABC matrix with ode45
                elseif input.AlgorithmType{2}==2 % ABC matrix with ode45
                    if input.reduced{i}==1
                        err('I still have to do all the projections on the eigenvector matrix')
                    else
                        n=length(d0);
                        A = [zeros(n,n), eye(n);-Step{step_counter}.M22\Step{step_counter}.K22, zeros(n,n)];

                        if input.discreteVScontinuous==1 % discrete force matrices

                            err('I do not think it is worth to spend time on a staet vector augmented solver for forces given as fixxed arays in time and do not take advantage of ode45 otherwise just use galpha and do not ask too much')
                        else % discreteVScontinuous==2 ... for continuously varing force vector with time
                            state_vec0=[d0;v0];
                            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                            % the problematic part %%%%%%%%%% I can write a
                            % generalized alpha solver but I am struggling
                            % with running embedded matlab functions
                            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                            f=@(state_vec,time)eqfunODE(n,A,Step{step_counter}.M22,Step{step_counter}.K22,state_vec,time,TotalGlobForceVecReduced{i});

                            [t,Y001] = ode45(@eqfunODE,[input.time_check_point{i},input.time_check_point{i+1}],state_vec0);
                            [t,Y002] = ode45(f,[input.time_check_point{i},input.time_check_point{i+1}],state_vec0); %differnet way of writing just to check if it is the same and it is 
                        end

                    end
                else
                end
        elseif input.StepSolverType{i}==3 % I might put a non linear static solver with geometrical stiff
        else
        end
        
        % use for continuity over loadsteps
%         t0=tc;
%         d0 = dM(:,end);
%         v0 = vM(:,end);
%         a0 = aM(:,end);
%         dMtot=cat(2,dMtot,dM);
%         vMtot=cat(2,vMtot,vM);
%         aMtot=cat(2,aMtot,aM);
    end
    % this uually are the output of this function
%     sol.dMtot=dMtot;
%     sol.vMtot=vMtot;
%     sol.aMtot=aMtot;
%     sol.Step=Step;
%     sol.timeline=timeline;


function [dydt]=eqfunODE(n,A,M22,K22,state_vec,time,TotalGlobForceVecReduced)
    y1=zeros(n,1);
    y2=zeros(n,1);
    dydt_1=zeros(n,1);
    dydt_2=zeros(n,1);
    dydt=zeros(2*n,1);
    y1=state_vec(1:n);
    y2=state_vec(n+1:2*n);
    dydt_1=y2;
    dydt_2=-M22\K22*y1+M22\TotalGlobForceVecReduced(time);
    dydt = [dydt_1; dydt_2];
    
    % other approach not working
%     Fexternal=[zeros(n,1);M22\TotalGlobForceVecReduced(time)];
%     
%     
%     
%     equdif=A*state_vec(1:2*n,1)+Fexternal(time);
end

Thank you in advice!!

Answer 1

[Ameer Hamza]

This call to ode45 is not correct:

[t,Y001] = ode45(@eqfunODE,[input.time_check_point{i},input.time_check_point{i+1}],state_vec0);

ode45 accepts function handle with two inputs. eqfunODE has 7 inputs. You can remove this line. The following lines are correct with a minor mistake

f=@(state_vec,time)eqfunODE(n,A,Step{step_counter}.M22,Step{step_counter}.K22,state_vec,time,TotalGlobForceVecReduced{i});
[t,Y002] = ode45(f,[input.time_check_point{i},input.time_check_point{i+1}],state_vec0); %differnet way of writing just to check if it is the same and it is 

You inverted the order of time and state_vec. The correct syntax should be

f= @(time,state_vec) eqfunODE(...

Now at least the syntax is correct, but I cannot run your code. Although you have attached, the function handles, e.g., TotalGlobForceVecReduced but they don't contain the definition of the external function they called, e.g., there is a call to cellfunc. You can try to run this and see if there is an error, you can show the error message.

[Francescogiuseppe Morabito]

eheheh.. how much one can feel ashamed when such things happens?

Anyway, thank you for the solutionand thank you for the advice on what is accepted as input. I am still a bit wandering though why the order is so important when the way it was defined in the function was state_vec,time. Does Ode45 have to get them in such order?

By the way now it seems that it is working correcty with no errors coming out but I still have to check if the solution is anywher correct, but I think I can nailed it. Thank you again and sorry for the function handle, I thought It was added to the script.

Just one last question.. I have definded the odefunction to give in output a column vect, is it normal that it is giving the output in strings? Sorry but this was my very first time with ode45.

thank you again and cheers

[Ameer Hamza]

No problem! Sometimes the syntax can be confusing when using a function first time.

Yes, the order is important for ode45. Because it treats your ODE function as a black box, it wants to input time and state vector and get the derivative at those states as output. It does not know the internal formula of how you calculate the derivative. So it must always give inputs to ODE function and get one output in a specific order.

What do you mean that it is giving output as a string? I guess output of ODE function should be numeric.

[Francescogiuseppe Morabito]

sorry I meant raw vector

[Ameer Hamza]

I am not sure why that would happen. From the definition of eqfunODE, output seems like a column vector.

[Francescogiuseppe Morabito]

mmm.. I will take a look at that. For now my state vector is 120x1 and the solution Y2 is 41x120 while t is 41 element array. Thank you for everything.

[Ameer Hamza]

Oh! I misunderstood in my previous comment. The ode45 returns the states as rows in the output matrix. The solution is correct. The ode45 retuned solution at 41 time values. Each row in Y2 is the solution at the time in the corresponding row in 't'.

[Francescogiuseppe Morabito]

ah okok I was guessing that the 41 size was reflecting the time steps. So everything is fine. Thank you!

[Ameer Hamza]

I am glad to be of help.