
Not getting six coupled ODEs solution using bvp4c in MATLAB, Please help me
13 views (last 30 days)
Show older comments
AISHWARYA KASARLA
on 9 Jun 2020
Commented: AISHWARYA KASARLA
on 10 Jun 2020
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))];
0 Comments
Accepted Answer
Stephan
on 10 Jun 2020
Edited: Stephan
on 10 Jun 2020
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:

More Answers (0)
See Also
Categories
Find more on Function Creation in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!