How to check the convergence and accuracy of the Problem using BVP4c
11 views (last 30 days)
Show older comments
Hi.
I want to check the convergence of the following BVP.
Problem_1()
function Problem_1
global alpha beta Pr Nb Nt
%eq1
alpha=1.2; beta=0.3;
%eq2
Pr=1.7; Nb=0.1; Nt=0.1;
solinit = bvpinit(linspace(0,6),[0 1 0 1 0]);
sol= bvp4c(@shootode,@shootbc,solinit);
eta = sol.x;
f = sol.y;
figure(1)
plot(eta,f(2,:));
hold on
end
function dydx = shootode(eta,f);
global alpha beta Pr Nb Nt
dydx = [f(2)
f(3)
(-((4/3)*alpha*beta*(1+2*eta)*f(3)^3 )-2*alpha*f(3)+f(2)^2-f(1)*f(3))
f(5)
(-2*f(5)-Pr*f(1)*f(5)-Pr*(1+2*eta)*(Nb*f(5)+Nt*f(5)^2 ))
];
end
function res = shootbc(fa,fb)
res = [ fa(1)-0
fa(2)-1
fa(4)-1
fb(2)-0
fb(4)-0
];
end
Thanks in advance.
2 Comments
Torsten
on 19 Nov 2023
What exactly do you want to do ? Check whether the solution given by bvp4c is correct ? Check whether you use enough grid points to meet a desired precision for the solution ?
Answers (1)
Torsten
on 19 Nov 2023
Edited: Torsten
on 19 Nov 2023
Inf = 3 seems to be sufficient:
global alpha beta Pr Nb Nt
%eq1
alpha=1.2; beta=0.3;
%eq2
Pr=1.7; Nb=0.1; Nt=0.1;
hold on
for i = 0:2
solinit = bvpinit(linspace(0,3+i,(3+i)*100+1),[0 1 0 1 0]);
sol= bvp4c(@shootode,@shootbc,solinit);
eta = sol.x;
f = sol.y;
plot(eta,f(2,:));
end
hold off
function dydx = shootode(eta,f);
global alpha beta Pr Nb Nt
dydx = [f(2)
f(3)
(-((4/3)*alpha*beta*(1+2*eta)*f(3)^3 )-2*alpha*f(3)+f(2)^2-f(1)*f(3))
f(5)
(-2*f(5)-Pr*f(1)*f(5)-Pr*(1+2*eta)*(Nb*f(5)+Nt*f(5)^2 ))
];
end
function res = shootbc(fa,fb)
res = [ fa(1)-0
fa(2)-1
fa(4)-1
fb(2)-0
fb(4)-0
];
end
1 Comment
Torsten
on 19 Nov 2023
Or maybe like this:
global alpha beta Pr Nb Nt
%eq1
alpha=1.2; beta=0.3;
%eq2
Pr=1.7; Nb=0.1; Nt=0.1;
for i = 0:3
solinit = bvpinit(linspace(0,3+i),[0 1 0 1 0]);
sol= bvp4c(@shootode,@shootbc,solinit);
x = linspace(0,3,100);
f(i+1,:,:) = deval(x,sol);
end
hold on
for i = 1:3
plot(x,squeeze(abs(f(i+1,2,:)-f(i,2,:))))
end
hold off
function dydx = shootode(eta,f);
global alpha beta Pr Nb Nt
dydx = [f(2)
f(3)
(-((4/3)*alpha*beta*(1+2*eta)*f(3)^3 )-2*alpha*f(3)+f(2)^2-f(1)*f(3))
f(5)
(-2*f(5)-Pr*f(1)*f(5)-Pr*(1+2*eta)*(Nb*f(5)+Nt*f(5)^2 ))
];
end
function res = shootbc(fa,fb)
res = [ fa(1)-0
fa(2)-1
fa(4)-1
fb(2)-0
fb(4)-0
];
end
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!