Efficiently Repetitive Solving of ODE

Hi all,

I'm working on a problem where I have to solve a nonlinear 2nd order ODE with initial conditions. I have to solve this ode for a set of constants (a, b, c) that coorespond to unique initial values and all are passed into my function as 1xN arrays. The constants don't change all that much, nor to the initial conditions, but will have slightly different solutions to the ode. The ode has no known symbolic solution so I feel I am stuck with numerically solving it N times.

At the moment my code works, but it is slow. Is there a more efficient way to approach this problem? Somehow leverage that my 1xN arrays do not change much and the solutions will be similar?

EDIT 11/22: The overall intent is to vary the parameters/condititions and compare the solutions. Physically, the ODE is a function of space only, but the parameters vary with time.

function value=some_function1(a,b,c,slope,x0)
%a, b, c, slope, x0 -> size of 1xN%

value=nan(size(a));
ODEFUN=@(x,Y,a,b,c) [Y(2);(a-b./(1-x).^3-c./(1-x).^2-Y(2)./x./sqrt(1+Y(2).^2)).*-(1+Y(2).^2).^(3/2)];
for ii=1:length(a)
    sol=ode113(@(x,y) ODEFUN(x,y,a(ii),b(ii),c(ii)),[x0(ii) 0],[0 slope(ii)]);
    %some important output value is calculated%
    value(ii)=some_function2(sol.x,sol.y(1,:))
end

end

Thank you!

[Star Strider]

I’m not certain that I competely understand the problem.

Is the intent to fit parameters in the ODE such that the integrated ODE is the objective function for the parameter estimation optimisation, or simply to vary the intital conditions (or other paramerers) in the ODE integration to do a number of different simulations and then plot (or otherwise compare) ther results?

.

[Brandon Murray]

Thanks - I'll make it a bit clearer.

The intent is the latter. The actual solution of the ODE for each set of parameters/conditions is used to find the arc length of the soluton (and then compare that).

Physically, the ODE is a function of space only, but the parameters vary with time.

[James Tursa]

Do you have a small time span where linearization might give you a faster means of finding approximate neighboring solutions?

Or maybe just feed all your ode problems to a parfor loop.

[Star Strider]

If this is a system of partial differential equations (functions of time and space), there are functions such as pdepe to integrate them and the Partial Differential Equation Toolbox devoted to solving them.

.

[Brandon Murray]

It appears parfor gives a minor increase in speed, I may go that route because it's an easy change. Linearization is also possible - Thanks for the suggestion. I'll give that a try and see If a little work there gives progress

As far as it being a PDE goes, I don't believe it is possible to formulate that way, I cannot easily specify the derivative of a,b,or c with time - it would involve a derivative of experimental measurements and that would be unreasonably noisy, and sometimes the derivative of a curve fit or something might not physically represent the data at hand

Thank you both for your suggestions

Answer 1

[Walter Roberson]

If you have two odes with the same time span, then

[t1, y1] = ode(first_ode, tspan, first_conditions);
[t2, y2] = ode(second_ode, tspan, second_conditions);

is, at least in theory, equivalent to

[t3, y3] = ode(@(t,y) [first_ode(t,y(1:length(first_conditions))); second_ode(t,y(length(first_conditions)+1:end))], tspan, [first_conditions, second_conditions])

And this implies that when first_ode and second_ode are the same code, that you can slightly rewrite the ode:

function dy = joint_ode(t, yM)
    y = reshape(yM, number_of_parameters, []);
    dy = zeros(size(y));
    
    %now calculate using ROW indexing, such as
    dy(1,:) = y(2,:).^2 - y(3,:);
    dy(2:end-1,:) = y(3:end,:);
    dy(end,:) = something;
    
    %now reshape to vector
    dy = dy(:);
end

[Brandon Murray]

Thank you, however in the case the tspan is changing, tspan=[x0(ii) 0], so I don't think this would work as-is but I will consider if I can reformulate the problem to obtain this form