Boundary value problem-bvp4c
2 次查看(过去 30 天)
显示 更早的评论
sepideh
2023-7-3
Hello
I am working on a paper by Di Federico et al 2012, they solved this by Mathematica@. I need to reproduce it in MATLAB and extend it for radial geometry.
Here is the problem :
I am trying to solve an ODE ( has been derived from an PDE by scaling and similarity solution process) which I have two values of it at the end of the interval [0,1] .
I know that as , the function tends to zero .
Also I know that as the derivative equals to .
as "bvp4c" needs the values in two different point, how should I write the residual function "bcfun" for this problem as I do not have two points ?
instead I have information about one end point and the slope near that point
This is the system of first order ODE :
also there is singularity near this point which they handle it by studying the problem in the neighborhood
the solution should have an answer like this :
采纳的回答
Torsten
2023-7-3
Use ode45 and integrate backwards:
22 个评论
sepideh
2023-7-4
Thank you
I saw the link. but there is a problem
in the example to sent to me the end points are known (y0=1) (y0=2.7)
in my case I only know that in y(1)=0 and I dont have any information about the starting point
I just know the slop which is y'= -2/3
Torsten
2023-7-4
Of course - since you divide by Y1 - you cannot set Y1 = 0 at x = 1 exactly. Use some small number instead of 0, e.g. 1e-8.
sepideh
2023-7-6
Dear Toesten
I wrote 2 codes, non of them work.
Would you please take a look at them and kindly correct them?
sepideh
2023-7-6
eps= 1e-10;
xmesh= linspace(0,1-eps,10);
solinit= bvpinit(xmesh,[1,0]);
sol=bvp4c(@bvpfun,@bcfun,solinit);
plot(sol.x,sol.y,'o')
function dydx=bvpfun(x,y)
dydx= [y(2),-(3*y(2)^2+x*y(2)+y(1))/(3*y(1))];
end
function res=bcfun(~,yb)
res=[yb(2)+2/3
yb(1)-(2/3*eps)];
end
Torsten
2023-7-6
xspan = [1 0];
y0 = [1e-8 -2/3];
fun = @(x,y)[y(2);-(3*y(2)^2+x*y(2)+y(1))/(3*y(1))];
[X,Y] = ode45(fun,xspan,y0);
plot(X,Y)
sepideh
2023-7-6
编辑:sepideh
2023-7-6
xspan = [1 0];
y0 = [1e-8 -1/3];
fun = @(x,y)[y(2);-(3*y(2)^2+x*y(2)+y(1))/(3*y(1))];
[X,Y] = ode45(fun,xspan,y0);
plot(X,Y(:,1))
hold on
syms F(Y) %analytical solution
F(Y)= 1/6*(1-Y^2);
fplot(F,[0 1],'o')
axis([0 1 0 1])
hold off
thank you very much
now the analytical and numerical results are perfectly match. you helped me a lot
sepideh
2023-7-6
编辑:sepideh
2023-7-6
I am working on Gravity Currents and the ensuing ODE has a analytical solution for the shape function as : (see Huppert 1982, Di Federico 2012)
this is the case for sudden or instantaneous release (
for the ODE can not solved analytically
here I made a mistake, actually the slope for this situation is -1/3 not -2/3 because now I am trying to solve the ODE for instantaneous release. the code that you gave me was perfectly match with the analytical solution (I correctted the slope) in this situation
when (constant injection) the ODE becomes :
with the sloe tends to -2/3
this can not solved analytically.
I am going to test the code for this equation now.
sepideh
2023-7-6
编辑:sepideh
2023-7-6
xspan = [1 0];
y0 = [1e-8 -1/3];
fun = @(x,y)[y(2);-(3*y(2)^2+x*y(2)+y(1))/(3*y(1))];
[X,Y] = ode45(fun,xspan,y0);
plot(X,Y(:,1))
hold on
syms F(Y)
F(Y)= 1/6*(1-Y^2);
fplot(F,[0 1],'o')
axis([0 1 0 1])
hold on
xspan = [1 0];
y0 = [1e-8 -2/3];
fun = @(x,y)[y(2);-(3*y(2)^2+2*x*y(2)-y(1))/(3*y(1))];
[X,Y] = ode45(fun,xspan,y0);
plot(X,Y(:,1),'--')
hold off
as you can see the constant injection affect the height of shape function (match with our governing equation) , this is what I expected based on literature.
thank you
sepideh
2023-9-9
Dear Torsten
I am trying to use the code for radial geometry but the numerical solution doesnot match the analytical one.
Would you please take a look at this and help me find out what is the problem? somehow in the mid interval its deviated
the ODE reads:
analytical solution is : (dots)
I used the base that you sent to me for backward method :
% plot analytical solution.
syms F(Y)
F(Y)= 1/8*(1-Y^2);
fplot(F,[0 1],'o')
axis([0 1 0 1])
hold on
% numerical solution for alpha=0
xspan = [1 0];
y0 = [1e-8 -1/4];
fun = @(x,y)[y(2);-(4*y(2)^2+(x)*y(2)+2*y(1))/(4*y(1))];
[X,Y] = ode45(fun,xspan,y0);
plot(X,Y(:,1))
ylabel('\Phi(\zeta)')
xlabel('\zeta')
hold off
Torsten
2023-9-9
syms x psi(x) f(x)
f = 1/8*(1-x^2);
ode = 0.5*psi+0.25*x*diff(psi,x)+diff(psi*diff(psi,x));
res = simplify(subs(ode,psi,f))
res(x) =
sepideh
2023-9-9
Thanks for your answer
I am trying to reproduce the results of an article, in that article the results are totally match.
sorry, I do not understand your solution, you want to show that there is a residual?
Torsten
2023-9-9
Yes. The analytical expression 1/8*(1-x^2) is not a solution of your differential equation.
sepideh
2023-9-11
yes you are right. I made a mistake in my modeling.
the correct ODE is:
now there is no residual.
but still there is a problem in the pragh.
xspan = [1 0];
y0 = [1e-10 -1/4];
fun = @(x,y)[y(2);-(2*x*y(1)+x^2*y(2)+4*x*y(2)^2+4*y(1)*y(2))/(4*x*y(1))];
[X,Y] = ode45(fun,xspan,y0);
plot(X,Y(:,1))
Torsten
2023-9-11
Yo divide by x and y(1). Thus you have to start with something small and different from 0 for y(1) in the y0-vector and you must integrate up to something small and different from 0 for xspan(2) in the xspan-vector.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)