Not getting six coupled ODEs solution using bvp4c in MATLAB, Please help me

1 次查看(过去 30 天)
I am solving six coupled first order ODEs using bvp4c in MATLAB, but it is not working...here is my code
function finalprogramme
thetaa = 0; thetab = pi/2 ;
solinit = bvpinit(linspace(thetaa,thetab,1000),[0.015 0 0 1 0 0]);
sol = bvp4c(@sixodes,@bc,solinit);
thetaint = linspace(thetaa,thetab,1000);
Sthetaint = deval(sol,thetaint);
plot(thetaint,Sthetaint([1 3],:));
function res = bc(ya,yb)
L = 70; u_star = 769.23e2; F=0.5*L;
res = [ ya(2); ya(3); ya(5); yb(2); yb(3); yb(5)-(-F/(2*L*u_star*H))];
function dydtheta = sixodes(theta,y)
R = 10; neu_star=0.3; alpha=(neu_star)/(1-2*neu_star); K=(12*(R^2)*(1-neu_star))/(H^2);
dydtheta = [y(2); ((5*K*alpha+4*alpha+3*K+3)/(1+K*(alpha+1)))*y(4)+((2*(1+2*alpha)*(K+1))/(1+K*(alpha+1)))*y(1)-R*((1-K*alpha-K)/(1+K*(alpha+1)))*y(6); y(4); ((-1)/(2*(2*alpha+1)))*((3*alpha+2)*y(2)-alpha*y(3)-alpha*R*y(5)); y(6); ((K+2)/(2))*(((-1)/(Ri))*y(2)+((1)/(R))*y(3)+y(5))];

采纳的回答

Stephan
Stephan 2020-6-10
编辑:Stephan 2020-6-10
I used nested functions - for 2 of the constants you did not provide values:
% Call the outer function
finalprogramme
%% Outer function
function finalprogramme
% You missed to tell us values for these both constants - used 1 for them
H = 1;
Ri = 1;
% All other Constants
thetaa = 0;
thetab = pi/2;
R = 10;
neu_star=0.3;
L = 70;
u_star = 769.23e2;
F=0.5*L;
alpha=(neu_star)/(1-2*neu_star);
K=(12*(R^2)*(1-neu_star))/(H^2);
% Solve the problem
solinit = bvpinit(linspace(thetaa,thetab,1000),[0.015 0 0 1 0 0]);
sol = bvp4c(@sixodes,@bc,solinit);
thetaint = linspace(thetaa,thetab,1000);
Sthetaint = deval(sol,thetaint);
plot(thetaint,Sthetaint([1 3],:));
%% Inner functions
function res = bc(ya,yb)
res = [ ya(2); ya(3); ya(5); yb(2); yb(3); yb(5)-(-F/(2*L*u_star*H))];
end
function dydtheta = sixodes(~,y)
dydtheta = [y(2); ((5*K*alpha+4*alpha+3*K+3)/(1+K*(alpha+1)))*y(4)+((2*(1+2*alpha)*(K+1))/(1+K*(alpha+1)))*y(1)-R*((1-K*alpha-K)/(1+K*(alpha+1)))*y(6); y(4); ((-1)/(2*(2*alpha+1)))*((3*alpha+2)*y(2)-alpha*y(3)-alpha*R*y(5)); y(6); ((K+2)/(2))*(((-1)/(Ri))*y(2)+((1)/(R))*y(3)+y(5))];
end
end
Result (using my constants for H & Ri) is:

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Solver Outputs and Iterative Display 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by