MATLAB Answers

Error in ODE45, must return a column vector

7 ビュー (過去 30 日間)
Brad Scott
Brad Scott 2021 年 2 月 23 日
コメント済み: Adam Danz 2021 年 2 月 26 日
%%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 件のコメント
Adam Danz
Adam Danz 2021 年 2 月 26 日
@Brad Scott, in response to your flag, it would be better to improve the question by editing it rather than reposting it.

サインインしてコメントする。

回答 (1 件)

James Tursa
James Tursa 2021 年 2 月 23 日
編集済み: James Tursa 2021 年 2 月 23 日
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 件のコメント
Walter Roberson
Walter Roberson 2021 年 2 月 26 日
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(:);

サインインしてコメントする。

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by