How to check the convergence and accuracy of the Problem using BVP4c
12 次查看(过去 30 天)
显示 更早的评论
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 个评论
Torsten
2023-11-19
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 ?
回答(1 个)
Torsten
2023-11-19
编辑:Torsten
2023-11-19
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 个评论
Torsten
2023-11-19
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
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Boundary Value Problems 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!