odextend initial guess error

1 次查看（过去 30 天）
구구
구구 2023-10-14
评论： 구구 2023-10-16
I am using ode15s to analyze the steady state of the system from the initial state (time t0) and then using odextend to analyze the external force acting on the system from time t1 to t2.
The problem is that depending on the frequency of the external force, when I use the odextend function, I get the error "Need a better guess y0 for consistent initial conditions." error pops up.
If I do a transient analysis from time t0 to t2 without using odextend, I have no problem, so what is the problem?
I have a lot of interpretation cases and I want to use odextend to reduce the overall interpretation time.
Thank you
  1 个评论
Torsten
Torsten 2023-10-14
I don't understand what you alternatively do instead of the transient analysis from time t0 to t2 without odextend.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Sam Chak
Sam Chak 2023-10-15
Ran in:
Hi @구구
I believe the issue may be related to inconsistent initial values. Under normal circumstances, the subsequent extended solution in odextend() should continue integrating from sol.x(:,end), unless a new set of initial values, x0new, is specified, which may be causing internal dynamic conflicts in the new odefcn2().
% generate solution from t0 to t1
tspan = [0 10];
x0 = [1 0]; % initial values
sol = ode45(@odefcn1, tspan, x0);
t = linspace(0, 10, 1001);
x = deval(sol, t);
plot(t, x(1,:)), hold on
% extend solution from t1 to t2 by injecting an external force
solext = odextend(sol, @odefcn2, 20);
t = linspace(10, 20, 1001);
x = deval(solext, t);
plot(t, x(1,:), 'r'), hold off
grid on, xlabel('t'), ylabel('x(t)')
title('Sinusoidal force is applied from 10 s to 20 s')
%% Mass-Spring-Damper, without force
function dxdt = odefcn1(t, x)
dxdt = zeros(2, 1);
force = 0;
dxdt(1) = x(2);
dxdt(2) = force - 2*x(2) - x(1);
end
%% Mass-Spring-Damper, with force
function dxdt = odefcn2(t, x)
dxdt = zeros(2, 1);
amp = 0.5;
freq = 2*pi/10;
xref = amp*sin(freq*t);
Dxref = freq*amp*cos(freq*t);
D2xref = - freq*freq*amp*sin(freq*t);
force = D2xref + 2*Dxref + xref;
dxdt(1) = x(2);
dxdt(2) = force - 2*x(2) - x(1);
end
  2 个评论
구구
구구 2023-10-16
thank you very much.
I found that I made a small mistake when applying the external force, so the y0 error occured.
Sam Chak
Sam Chak 2023-10-16
Ran in:
You are welcome, @구구. In fact, since you mentioned that you want to test the responses of the system by injecting a range of frequencies of the external force, then @Torsten's for-loop approach is very efficient. For example:
freq = linspace(1, 10, 7)
freq = 1×7
1.0000 2.5000 4.0000 5.5000 7.0000 8.5000 10.0000
numFreq = numel(freq)
numFreq = 7
% perform iterations 7 times
for j = 1:numFreq
% insert the ode45 solver here
xref = amp*sin(freq(j)*t); % <-- find the j-th position in the array of 'freq'
end

请先登录，再进行评论。

更多回答（1 个）

Torsten
Torsten 2023-10-15
移动：Torsten 2023-10-15
Ran in:
  1. Compute the solution up to t1 using ode15s.
  2. Make a loop over the different F(t)'s you want to prescribe and restart all the continuing integrations using the normal call to ode15s from the solution obtained in step 1.
No need to use odextend.
a = 1;
fun = @(t,y) a*y;
t0 = 0;
t1 = 1;
t2 = 2;
y0 = 1;
[Tstart,Ystart] = ode15s(fun,[t0 t1],y0);
anew = [2 3 4];
hold on
for i = 1:numel(anew)
fun = @(t,y) anew(i)*y;
[Tcont,Ycont] = ode45(fun,[t1 t2],Ystart(end,:));
plot([Tstart(1:end-1);Tcont],[Ystart(1:end-1,:);Ycont]);
end
hold off
grid on

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by