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!!
采纳的回答
Ameer Hamza
2020-5-3
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.
8 个评论
Francescogiuseppe Morabito
2020-5-3
Ameer Hamza
2020-5-3
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
2020-5-3
Ameer Hamza
2020-5-3
I am not sure why that would happen. From the definition of eqfunODE, output seems like a column vector.
Francescogiuseppe Morabito
2020-5-3
Ameer Hamza
2020-5-3
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
2020-5-3
Ameer Hamza
2020-5-3
I am glad to be of help.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
