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...
clcclear allclose allload('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);
endThank you in advice!!
Best Answer