How to compare two separate outputs of the ode45 function?

So, I am trying to perform a parametric analysis using two separate differential equations. I am doing it by taking the xvalues at times 0 to 2 from one ode and subtracting from the other ode that has two independent variables.

This is the function code for one ode that has a changing time variable and changing ceq variable.

function xdot=linear(t,x)
ceq=1:.1:100;
k=15000;
m=10;
xdot(1,1)=x(2);
xdot(2,1)=-(ceq/m)*x(2)-(k/m)*x(1);

This is the function code for the other ode that has just time as the independent variable.

function xdot=funnonlinear(t,x);
c=38.11;
k=15000;
m=10;
u=.1;
g=9.81;
xdot(1,1)=x(2);
xdot(2,1)= -(c/m)*x(2)-(u*g)*sign(x(2))-(k/m)*x(1);

Here is the main code where I am trying to estimate ceq by comparing the average error of each output for each time for the nonlinear ode and output for each time of each different ceq values. I am also wanting the program to display ceq value of the smallest error output.

clear all; close all; clc;
t0=0;
tf=2;
tspan=[t0 tf];
x0=[0 1];
[ta xa]=ode45('funnonlinear',tspan,x0);
[ta xb] =ode45('linear',tspan,x0);
error=mean(abs(xb-xa));
for min(error)
    disp(ceq)
end

I am also getting this error:

Unable to perform assignment because the size of the left side is 1-by-1 and the size of the right side

is 1-by-991.

Error in linear (line 6)

xdot(2,1)=-(ceq/m)*x(2)-(k/m)*x(1);

Error in odearguments (line 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in CaseStudyP4 (line 7)

[ta xb] =ode45('linear',tspan,x0);

[Walter Roberson]

You have ceq = 1:.1:100 . That is 991 elements long. You use all of ceq in defining xdot(2,1) so you are trying to store 991 elements into a slot reserved for 1 element.

You will need to define function linear as expecting 992 boundary conditions instead of 2 boundary conditions, and you will need to use

xdot(2:992,1) = -(ceq/m)  etc

[Shawn Alcorn]

I fixed that part just now, but now I'm getting the error of the matrix dimesions not agreeing. xa is 158885x2 and xb is 437 x 992.

[Walter Roberson]

Pass in a vector of times for tspan to return the samples at only those times.

However, I have no idea on the intended meaning of comparing 2 outputs to 992 outputs.

[Shawn Alcorn]

What it is suppose to do is for each ceq value subtract the output of the linear function at each time from the output of the nonlinear function at each time and average the difference until a minimum value of the difference is found then display the ceq value that gave the minimum difference.

Answer 1

[Walter Roberson]

Sounds like a boundary volume problem to me. But perhaps you have a reason to be more or less using a shooting method.

function xdot=linear(t, x, ceq)
k=15000;
m=10;
x1 = reshape(x(1:2:end), [], 1);
x2 = reshape(x(2:2:end), [], 1);
nc = length(ceq);
txdot = zeros(2, nc);
txdot(1, :) = x2;
txdot(2, :) = -(ceq(:)/m) .* x2 - (k/m) .* x1;
xdot = txdot(:);

To be called like,

ceq = 1:.1:100;
nc = length(ceq);
x0 = [0 1];
[ta, xa] = ode45(@(t,x) funnonlinear(t, x, ceq), tspan, repmat(x0, 1, nc));

The return in xa will have 2*length(ceq) columns, and those columns will alternate, with odd numbered columns corresponding to xdot(1) and even numbered columns corresponding to xdot(2), with one pair of columns for each ceq value. Effectively instead of running a loop processing each ceq value separately, all of the ceq values are processed at the same time.

After that you would do something like

xa3 = reshape(xa, size(xa,1), 2, []);
err = mean(sum(abs(bsxfun(@subtract, xb, xa)),2),1);
[min_err, min_idx] = min(err);
best_ceq = ceq(min_idx);