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!