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