ODE45 related

Hi all I have a ODE of the form:

 y"+ay'+by = 0 ------ (1)

 For a given time span of [0 t] and initial values, using ODE45 (or any solver) one can find the numerical solution to Eq. (1).

Ex: y = ode45('eqn',[0 t],[0 1]) ;

y is a solution to Eq. (1). It contains y in the first column, y' in the second column. I want to find the derivative of y w.r.t to the coefficients of Eq. (1). i.e I want to find dy/da and dy/db. How to do this?

Thanks in advance

Sreenu

Answer 1

[Jan]

You can vary a by a sufficiently small number for each component and recalcute the trajectory. Then you get the derivative of the trajectory to the paremeter by the quotient of differences. This is called external numerical differentiation and it has the drawback, that "sufficiently small" is not well-defined.

For an internal numerical differentiation you have to vary the trajectory in each step of the solver - as far as I know this is not possible with the original ODE45. Expanding ODE45 is not really hard, but finding the best width for the local variation in each step remains a magic piece of science. Actually at least an approximation of the 2nd derivatives are required for valid variation widths - but this is very expensive.

There must be an integrator in Matlab, which calculates an internal numerical differentiation, because this is needed for the optimization. Perhaps another user knows how to manage this.

Answer 2

[AL]

Can I ask one doubt in this only? I am also doing something similar. Solving 4th degree of freedom vibration system using following formula:

while using this approch I am getting 0 as phase value for real solution and for imaginary solution i am getting negative phase value; however, for real solution values should be greater then zero. Can anyone help me with this?

my force and damping frequancy are a function of time.

this is the code:

% Define the mass, damping, and stiffness matrices

load C.mat

load K.mat

load M.mat

% Define the input force as a sinusoidal function of frequency w in Hz

f = linspace(0, 40, 100); % linearly spaced frequency vector in Hz

w = 2*pi*f; % angular frequency vector in rad/s

F = @(t) [sin(2*pi*f(1)*t); zeros(3, 1)];

% Define the initial conditions

x0 = zeros(4, 1);

v0 = zeros(4, 1);

y0 = [x0; v0];

ind = 100;

% Define the time span

tspan = [0 20];

% Preallocate arrays to store the magnitude and phase of the response

mag = zeros(size(w));

phase = zeros(size(w));

X = zeros(length(w), 4);

% Loop over the frequency vector and compute the response at each frequency

for i = 1:length(w)

% Define the input force at the current frequency

F = @(t) [sin(w(i)*t); zeros(3, 1)];

% Solve the differential equation to obtain the steady-state response

[~, y] = ode45(@(t, y) vibration_equation(t, y, M, C, K, F), tspan, y0);

x = y(end, 1:4)';

% Compute the magnitude, phase, and displacement of the response at the current frequency

mag(i) = abs(x(1));

phase(i) = angle(x(1)) * 180/pi;

X(i,:) = x';

mpp = max(mag(1,:))/30;

end

% find the 4 dominant peaks

[Ypk,Xpk,Wpk,Ppk] = findpeaks(mag(1,:),'NPeaks',4,'SortStr','descend');

Xpk = Xpk + min(size(mag)) - 1;

figure;

subplot(2, 1, 1);

%plot(f, mag);

plot(f,mag,f(Xpk),Ypk,'dr');

xlabel('Frequency (Hz)');

ylabel('Magnitude');

ylim([0 3.5/10000]);

title('Frequency Response Function');

grid on;

subplot(2, 1, 2);

plot(f, mod(phase+180,360)-180,f(Xpk),Ypk,'dr');

xlabel('Frequency (Hz)');

ylabel('Phase (deg)');

ylim([-180, 180]);

grid on;

clear C damp_ratio F i ind K K1_Th K2_Th M Natural_freq tspan mxidx

function dydt = vibration_equation(t, y, M, C, K, F)

% Compute the acceleration

x = y(1:4);

v = y(5:8);

a = M \ (F(t) - C*v - K*x);

% Compute the derivative of the state vector

dydt = [v; a];

end

Current results:

I want something like this: