How to compare two separate outputs of the ode45 function?

5 次查看（过去 30 天）
Shawn Alcorn
Shawn Alcorn 2019-3-9
回答： Walter Roberson 2019-3-9
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);
  4 个评论
显示 2更早的评论
Walter Roberson
Walter Roberson 2019-3-9
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
Shawn Alcorn 2019-3-9
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.

请先登录，再进行评论。

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

回答（1 个）

Walter Roberson
Walter Roberson 2019-3-9
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);

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by