odextend initial guess error

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

[Torsten]

I don't understand what you alternatively do instead of the transient analysis from time t0 to t2 without odextend.

Answer 1

[Sam Chak]

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

[구구]

thank you very much.

I found that I made a small mistake when applying the external force, so the y0 error occured.

[Sam Chak]

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

Answer 2

[Torsten]

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

[구구]

thank you. I'm gonna try yhis