How to check the convergence and accuracy of the Problem using BVP4c

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.

[Torsten]

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 ?

[Syed Sohaib Zafar]

I want to do an analysis for the different grid points and compare them. To show that what's an appropriate selection for the infinity. Like in the above code I choose '6'.

Answer 1

[Torsten]

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

[Torsten]

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