Error in ODE45, must return a column vector
2 views (last 30 days)
Show older comments
%%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
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.
Answers (1)
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
on 26 Feb 2021
[t, y] = Cbar;
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
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(:);
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!