MATLAB Answers

Keeping values of function when solving ODE

조회 수: 1(최근 30일)
GCats
GCats 2021년 5월 6일
댓글: Bjorn Gustavsson 2021년 5월 6일
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?
Thank you in advance!

채택된 답변

Bjorn Gustavsson
Bjorn Gustavsson 2021년 5월 6일
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
Bjorn Gustavsson
Bjorn Gustavsson 2021년 5월 6일
Good, then I interpreted your requirements "correct enough".

댓글을 달려면 로그인하십시오.

추가 답변(0개)

태그

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by