# Keeping values of function when solving ODE

1 view (last 30 days)
GCats on 6 May 2021
Hello everyone!
I have a function describing a SDOF state space with a sinusoidal input y in one script as follows:
function dz = system(t,z,option,K,M,fe,y0)
dz=zeros(2,1);
y = y0 * sin(2*pi*fe*t); %base excitation
dz(1)=z(2);
dz(2)=(K*(y-z(1))/M
end
while in a separate script, I solve the ODE for different values of fe in a loop as follows:
clc
clear all
close all
% Linear Parameters________________________________________________________
K=5000;% linear stiffness
M=1; % system's mass
%Input Parameters__________________________________________________________
y0 = 1; % amplification of base input
T = 10;
fe = [1 2 3 4];
n_freq = length(fe);
for kk=1:n_freq
[tn,z] = ode15s('systeme1DDL_Dahl_bal',[t(1) t(end)],K,M,y0,fe(kk));
end
How can I save the value of y (the base excitation input) for each iteration?

Bjorn Gustavsson on 6 May 2021
As your ODE is written above we have to ask: Why bother? y is deterministically given provided your input parameters and the time. It should be trivial to calculate it for the times in you tn output. In case you have a slightly more complicated driving function you could perhaps do something like this:
function dz = system(t,z,option,K,M,fe,y_fcn)
dz=zeros(2,1);
y = y_fcn(t); %base excitation
dz(1)=z(2);
dz(2)=(K*(y-z(1))/M;
end
This would let you send in "any" function to your ODE:
% Linear Parameters________________________________________________________
K=5000;% linear stiffness
M=1; % system's mass
%Input Parameters__________________________________________________________
y0 = 1; % amplification of base input
T = 10;
fe = [1 2 3 4];
n_freq = length(fe);
z0 = 0; % you'll have to set this
for kk=1:n_freq
f_current = @(t) y0*sin(2*pi*fe(kk)*t);
[tn,z] = ode15s(@(t,z)system(t,z,K,M,fe(kk),f_current),[t(1) t(end)],z0);
y_exc{kk} = f_current(tn);
end
HTH
##### 2 CommentsShowHide 1 older comment
Bjorn Gustavsson on 6 May 2021
Good, then I interpreted your requirements "correct enough".