Vicious Circle Problem: How to solve ODEs with variables that are based on the solutions of the ODEs?

Assuming I have the following system of ODEs:

dSdt = [S(2);
        L*S(1)+R*(S(1)-S(3));
        S(4);
        -g*(S(1)-S(3))];

and the tspan

tspan = [0 t_max];

I want to solve the odes with ode45().

[t,S] = ode45(odefcn,tspan,y0)

The Probelm is, that the factor L depends on the outcomes of S(1) and S(3)

For example L could be:

S1_tmax = S(end,1);
S3_tmax = S(end,3);

L = S1_tmax + S3_tmax;

where S_max(1) and S_tmax(3) are the datapoints of S(1) and S(3) at the last possible points in time of the solution at t_max.

Is that even achievable? How can I break out of this vicious circle of cause and condition?

Answer 1

[Bjorn Gustavsson]

Perhaps something like this:

odeFCN = @(t,S,L,R) [S(2);
                   L*S(1)+R*(S(1)-S(3));
                   S(4);
                  -g*(S(1)-S(3))];
t_span = [t_max 0];
S_end0 = [9 8 7 5];
R = exp(pi);
S_0target = [1 2 3 4]; % You have some initial condition at t=0 I guess;
L = S_end0(1)+S_end0(3); % Whatever you have for this
[t,S] = ode45(@(t,S) odeFCN(t,S,L,R),t_span,S_end0);
while sum(abs(S(1,:)-S_0target)) > 1e-6
    S_end0 = % Some cunning update of the end boundary condition
    L = S_end0(1)+S_end0(3); % Whatever you have for this
    [t,S] = ode45(@(t,S) odeFCN(t,S,S_end0(1)+S_end0(3),R),t_span,S_end0);
end

HTH

[e_frog]

Thanks for the help so far!

Sadly your code stops working on this line:

Index exceeds the number of array elements (1).

Error in test (line 13)
    L = S_end0(1)+S_end0(3); % Whatever you have for this
    

Also I quite dont understand what you mean by this line:

    S_end0 = % Some cunning update of the end boundary condition

S_end is not known but has to be calculated. Thats the oint of my Problem. I need S_1end and S_3end to calculate L, but to get S_1 and S_3 the variable L has to be known.

[Bjorn Gustavsson]

You have to specify your problem well enough to be solvable. As of now you have told us that you have 4 coupled ODEs, where one coefficient will depend on the final solution. (This oviously breaks the causality of physics, and mathematically it doesn't make much sense as an integration with respect to time - what says that the coefficient is determined at an end-time of any given t_end, and not a bit later?) In order to make this a problem that's not trivially solvable (Until you give us an initial condition, we can just pick any random end-solution and calculate the coefficient L from that and integrate the ODEs backwards in time, and be done with it. That would give one solution of the coupled ODEs satisfying the request that L depends on the end-solution in the way you specify.

For this to be a solution with one (or a few) solutions you will have to specify a solution at t=0 that can be used as a target for a shooting-method, starting at t_end. Then the task is to search for values of S(t_end) that will give you a solution that hits S(t_0). Since you now start at t_end you will have values of S at that point, the task now has become "find S(t_end) that gives you S(t_0) with these ODEs", that is the shooting method - here applied in the reverse time-direction.

You have to do this task not me.

Answer 2

[Bjorn Gustavsson]

You might run the ODE-integration backwards, that would give you a problem where you could determine your unknown parameters from the initial conditions (at t_max). If you have initial conditions at t0, then it seems you have a problem that could be solved with the shooting method.

HTH

[e_frog]

Could you please help me figure out how to write the ode45() for my problem with the shooting method? I am struggeling to put in the inital conditions, tspan and the boundary conditions.

tspan = [0 t_max];
init = [0 2 0 2];
A = L - S(3);
b = L - S(1);
bndcnd = [t_max A t_max B];

[t,S] = ode45(@odefcn, ??????)

[e_frog]

Actually I learned from these file exchange comments, that the shooting method can not be used for a system of equations. Could you please help me write the matlab code for my problem?

[Bjorn Gustavsson]

Google and you will find. You haven't even specified any boundary conditions whatsoever yet...