Error in ODE45, must return a column vector

3 views (last 30 days)
Brad Scott
Brad Scott on 23 Feb 2021
Commented: Adam Danz on 26 Feb 2021
%%Can find the cbar solution using quadratic formula
%%Constants
D = (10^(-4))/(1-10^(-4));
gamma = 0.4;
M = 0.2;
z = linspace(0,1,1000);
mcQ = 10^5; %mathcal Q
n=2;
w=0.5;
%%Quadratic coefficients
a = 1;
b = D+gamma*z-M;
c = -D*M;
%%Quadratic formula, c is negative, so for real solution we just use + root
cbar = 0.5*(-b+sqrt(b.^2-4*a*c));
%Q from definition
Qbar = gamma*z + cbar - M;
%%Test solutions, avg error of -2.67*10^(-19) approx 0
test_c_Q=mean(cbar.*(D+Qbar) - D*M);
%%Find phi using root finder
for Qs=Qbar
root_est = (Qs/mcQ)^(1/n);
p(Qs==Qbar) = fzero(@(p) fphi(p,n,mcQ,Qs), root_est);
end
%%Test phi value, avg error of -6.46*10^(-17) approx 0
test_p=mean(p+mcQ*p.^(n).*(1-p).^2 - Qbar);
%%Setup to solve ODE
dQdp=1+n.*mcQ.*pbar.^(n-1).*(1-pbar).^2 - 2.*mcQ.*pbar.^2.*(1-pbar);
bdQdz=gamma./(1+D*M./(D+Qbar).^2);
bdcdz=-D*M./(D+Qbar).^2.*bdQdz;
ode = @(Qhat,chat) [1i.*w.*(X(1)./dQdp - X(2)) + X(2)./(D+pbar+cbar).*(1i.*w.*(D+pbar+cbar)-bdQdz)-X(1).*bdcdz, ...
X(2)./(D+Qbar+cbar).*(1i*w.*(D+pbar+cbar)-bdQdz)-X(1).*bdcdz];
ic = [bdQdz bdcdz+gamma];
[t,X] = ode45(ode,[0 1],ic);
function y = fphi(p,n,mcQ,Q)
y = p+mcQ*p^n*(1-p)^2-Q ;
end
Hey there!
I'm trying to solve the ODE above using the above IC's. Can anyone offer some help?
  1 Comment
Adam Danz
Adam Danz on 26 Feb 2021
@Brad Scott, in response to your flag, it would be better to improve the question by editing it rather than reposting it.

Sign in to comment.

Answers (1)

James Tursa
James Tursa on 23 Feb 2021
Edited: James Tursa on 23 Feb 2021
Just make your function handle return a column vector by using ; instead of , to separate the elements. E.g.,
ode = @(Qhat,X) [1i.*w.*(X(1)./dQdp - X(2)) + X(2)./(D+pbar+cbar).*(1i.*w.*(D+pbar+cbar)-bdQdz)-X(1).*bdcdz; ...
X(2)./(D+Qbar+cbar).*(1i*w.*(D+pbar+cbar)-bdQdz)-X(1).*bdcdz];
What is pbar?
  5 Comments
Walter Roberson
Walter Roberson on 26 Feb 2021
[t, y] = Cbar;
ans = 1×2
1 1
ans = 1×2
1 1000
ans = 1×2
1 1
ans = 1×2
1 1000
ans = 1×2
1 1000
ans = 1×2
1 1000
ans = 1×2
1 1000
ans = 1×2
2 1000
Error using odearguments (line 93)
CBAR/NESTED_ODEFUN must return a column vector.

Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in solution>Cbar (line 60)
[t,X] = ode45(@nested_odefun,[0 1],ic);
plot(t,y)
function [t, y] = Cbar
%%Can find the cbar solution using quadratic formula
%%Constants
D = (10^(-4))/(1-10^(-4));
gamma = 0.4;
M = 0.2;
z = linspace(0,1,1000);
mcQ = 10^5; %mathcal Q
n=2;
w=0.5;
%%Quadratic coefficients
a = 1;
b = D+gamma*z-M;
c = -D*M;
%%Quadratic formula, c is negative, so for real solution we just use + root
cbar = 0.5*(-b+sqrt(b.^2-4*a*c));
%Q from definition
Qbar = gamma*z + cbar - M;
%%Test solutions, avg error of -2.67*10^(-19) approx 0
test_c_Q=mean(cbar.*(D+Qbar) - D*M);
%%Find phi using root finder
for Qs=Qbar
root_est = (Qs/mcQ)^(1/n);
pbar(Qs==Qbar) = fzero(@(p) fphi(p,n,mcQ,Qs), root_est);
end
%%Test phi value, avg error of -6.46*10^(-17) approx 0
test_p=mean(pbar+mcQ*pbar.^(n).*(1-pbar).^2 - Qbar);
%%Setup to solve ODE
dQdp=1+n.*mcQ.*pbar.^(n-1).*(1-pbar).^2 - 2.*mcQ.*pbar.^2.*(1-pbar);
bdQdz=gamma./(1+D*M./(D+Qbar).^2);
bdcdz=-D*M./(D+Qbar).^2.*bdQdz;
function ode = nested_odefun(t, X)
size(w)
size(dQdp)
size(D)
size(pbar)
size(cbar)
size(bdQdz)
size(bdcdz)
ode(1,:) = 1i.*w.*(X(1)./dQdp - X(2)) + X(2)./(D+pbar+cbar).*(1i.*w.*(D+pbar+cbar)-bdQdz)-X(1).*bdcdz;
ode(2,:) = X(2)./(D+Qbar+cbar).*(1i*w.*(D+pbar+cbar)-bdQdz)-X(1).*bdcdz;
size(ode)
end
ic = [bdQdz(1) bdcdz(1)+gamma];
[t,X] = ode45(@nested_odefun,[0 1],ic);
function y = fphi(p,n,mcQ,Q)
y = p+mcQ*p^n*(1-p)^2-Q ;
end
end
Walter Roberson
Walter Roberson on 26 Feb 2021
My interpretation:
You are trying to solve a system of 1000 differential equations, but you only initialize with two boundary conditions.
If you passed in 2000 boundary conditions, interleaved, then at the end of the nested_odefun that I put in, put in
ode = ode(:);

Sign in to comment.

Tags

Community Treasure Hunt

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

Start Hunting!