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!