Error in ODE45, must return a column vector

%%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?

[Adam Danz]

@Brad Scott, in response to your flag, it would be better to improve the question by editing it rather than reposting it.

Answer 1

[James Tursa]

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?

[dpb]

But ode uses variables Qhat,chat as its dummy arguments that aren't used in the definition -- and does use X that aren't arguments???

[James Tursa]

Yes, I made an assumption that the 2nd argument should have been X. But that does raise the question you asked. It doesn't seem everything is defined/used properly.

[Brad Scott]

Hello there,

First of all thank you for responding. I had variable p changed to pbar, just forgot to do it in the code!

Basically I have two ODE's involving Qhat and chat, and I looked online at it said to replace those with X(1) and X(2) for the variables. I've tried replicating that in my code. Here's an updated version with another error.

%%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;

ode = @(t,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];

ic = [bdQdz(1) bdcdz(1)+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

[Walter Roberson]

[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]

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(:);