Info

此问题已关闭。 请重新打开它进行编辑或回答。

ODE45 related

2 次查看（过去 30 天）
KSSV
KSSV 2012-1-23
回答： AL 2023-3-21
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

此问题已关闭。

采纳的回答

Jan
Jan 2012-1-23
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.

更多回答（1 个）

AL
AL 2023-3-21
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:

此问题已关闭。

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by