Error in ODE45, must return a column vector
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
%%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 commentaire
Adam Danz
le 26 Fév 2021
@Brad Scott, in response to your flag, it would be better to improve the question by editing it rather than reposting it.
Réponses (1)
James Tursa
le 23 Fév 2021
Modifié(e) : James Tursa
le 23 Fév 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 commentaires
Walter Roberson
le 26 Fév 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
le 26 Fév 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(:);
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!