I'm trying too solve a set of coupled ODEs using BVP5C but I keep getting the same (wrong) results for the different cases. What am I doing wrong?
1 view (last 30 days)
Show older comments
I'm trying to solve a set of coupled ODEs using BVP4C but the results I get from Matlab isn't making sense. I keep getting the same value.
n = 0.4;
Bo = [0.01322917839
0.01082513952
0.004466835920
0.001258925411
0.0003548133892];
dy = zeros(5,1);
dy(1) = y(2);
dy(2) = ((2*n+1)*y(3)*y(2)*Bo^(2/(n+1))/(n+1))/0.2e1;
dy(3) = y(4);
dy(4) = y(5);
dy(5) = ((-y(4)^2*n+2*y(3)*y(5)*n+2*y(1)*n-y(4)^2+y(3)*y(5)+2*y(1))/y(5)^(n-1)/n/(n+1))/0.2e1;
Where,
y(1) = theta
y(2) = theta'
y(3) = f
y(4) = f'; and
y(5) = f''
Subject to the following Boundary Conditions,
3 initial conditions are given: eta=0, f(0)= f'(0)=0,theta(0)=1
The boundary conditions that needs to be satisfied are: f'(eta=inf)= 0 and theta(eta=inf)=0.
I have tried the following code but I get f''= 0 and theta'= -0.033333333333333 for all cases.
function sol= proj
clc;clf;clear;
global Bo n
% Bo=1;
% n=1;
pp=[0.4];
qq=[0.01322917839
0.01082513952
0.004466835920
0.001258925411
0.0003548133892];
%figure(1)
%plot(2,1);hold on
options=bvpset('stats','on','RelTol',1e-9);
m=linspace(0,30);
solinit= bvpinit(m,[1,0,0,0,0]);
for i=1:numel(pp)
n=pp(i);
for i=1:numel(qq);
Bo=qq(i)
sol= bvp5c(@projfun,@projbc,solinit,options);
format long
y1=deval(sol,0)
solinit= sol;
plot(sol.x,sol.y(1,:));hold on
end
end
end
function f= projfun(x,y)
global Bo n
f = zeros(5,1);
f = zeros(5,1);
f(1) = y(2);
f(2) = ((2*n+1)*y(3)*y(2)*Bo^(2/(n+1))/(n+1))/0.2e1;
f(3) = y(4);
f(4) = y(5);
f(5) = ((-y(4)^2*n+2*y(3)*y(5)*n+2*y(1)*n-y(4)^2+y(3)*y(5)+2*y(1))/y(5)^(n-1)/n/(n+1))/0.2e1;
end
function res= projbc(ya,yb)
res= [ya(1)-1; ya(3); ya(4); yb(1); yb(4)];
end
The code works perfectly for the following conditions when n =1 but gives the same wrong values for n other than 1.
n = 1.0;
Bo = [0.75
0.5
0.1
0.01
0.001
];
dy = zeros(5,1);
dy(1) = y(2);
dy(2) = ((2*n+1)*y(3)*y(2)*Bo^(2/(n+1))/(n+1))/0.2e1;
dy(3) = y(4);
dy(4) = y(5);
dy(5) = ((-y(4)^2*n+2*y(3)*y(5)*n+2*y(1)*n-y(4)^2+y(3)*y(5)+2*y(1))/y(5)^(n-1)/n/(n+1))/0.2e1;
What I'm I doing wrong?
10 Comments
Answers (0)
See Also
Categories
Find more on Boundary Value Problems 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!