Borrar filtros
Borrar filtros

how to write matlab bvp4c code for moving boundary surface problem?

9 visualizaciones (últimos 30 días)
hello, I am new in the bvp4c method of Matlab. I want to find multiple graphs for different value of 'r' for the following problem. Equation is f''' +ff''=0 Boundary condition is f(0)=0, f'(0)=1-r, f'(infnty)=r

Respuesta aceptada

Stephan
Stephan el 6 de Nov. de 2018
Editada: Stephan el 6 de Nov. de 2018
Hi,
i recommend to start simple and work step by step:
  1. Have a read in the documentation of bvp4c
  2. Try to implement your model with a single fixed value for r and get a valid solution
  3. Once you have this, use a for loop to vary r like needed for your purposes and save all the results in a Matrix in different columns or rows
  4. Plot the results by using the saved values
Problems by doing this? --> Come back and ask a specific question including your code so far.
Best regards
Stephan
  2 comentarios
Anil Gautam
Anil Gautam el 7 de Nov. de 2018
Editada: Torsten el 7 de Nov. de 2018
function cortell
global r
r=[0]
solinit = bvpinit(linspace(0,10,1000),[0 1- r r]);
options = bvpset('stats','on');
sol = bvp4c(@equation,@bvpbc,solinit,options);
eta = sol.x;
u = sol.y;
clf reset
plot(eta,u(3,:));
hold on;
axis([0 4 0 1.3]);
xlabel('\eta')
ylabel('f''(\eta)')
shg
% --------------------------------------------------------------------------
function dy=equation(~,y)
%dy=zeros(4,1);
dy=[y(2);
y(3);
-y(1)*y(3);
%---------------------------------------------------------------------
];
function res = bvpbc(y0,yinf)
global r
res = [y0(1)
y0(2)-1+ r
yinf(2)- r
];
if true
% code
end
this is code . how I can I get the different graph for different value of 'r' . In other word, for r=0.0, 0.25,0.75 how i can plot graph for f on y axis and eta on x axis.. please suggest.
Stephan
Stephan el 7 de Nov. de 2018
See the answer from Torsten...

Iniciar sesión para comentar.

Más respuestas (1)

Torsten
Torsten el 7 de Nov. de 2018
Editada: Torsten el 7 de Nov. de 2018
function main
r=[0 0.02 0.04]
for i=1:numel(r)
r_actual = r(i);
solinit = bvpinit(linspace(0,10,10000),[0 1-r_actual 0]);
options = bvpset('stats','on');
sol{i} = bvp4c(@equation,@(y0,yinf)bvpbc(y0,yinf,r_actual),solinit,options);
end
for i=1:numel(r)
eta = sol{i}.x;
u = sol{i}.y;
plot(eta,u(3,:));
hold on;
end
end
% --------------------------------------------------------------------------
function dy=equation(~,y)
%dy=zeros(4,1);
dy=[y(2);
y(3);
-y(1)*y(3);
%---------------------------------------------------------------------
];
end
function res = bvpbc(y0,yinf,r_actual)
res = [y0(1)
y0(2)-1+ r_actual
yinf(2)- r_actual
];
end
  3 comentarios
Torsten
Torsten el 7 de Nov. de 2018
Editada: Torsten el 7 de Nov. de 2018
Try this code - it uses the solution of the last step as initial guess for the next:
function main
global ix
r=[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0];
r_actual = r(1);
ir=1;
solinit = bvpinit(linspace(0,10,1000),@(x)guess(x,ir,[0 1-r_actual r_actual]));
options = bvpset('stats','on');
ix = 0;
sol{1} = bvp4c(@equation,@(y0,yinf)bvpbc(y0,yinf,r_actual),solinit,options);
for ir=2:numel(r)
r_actual = r(ir);
solinit = bvpinit(sol{ir-1}.x,@(x)guess(x,ir,sol{ir-1}.y));
options = bvpset('stats','on');
ix = 0;
sol{ir} = bvp4c(@equation,@(y0,yinf)bvpbc(y0,yinf,r_actual),solinit,options);
end
for ir=1:numel(r)
eta = sol{ir}.x;
u = sol{ir}.y;
plot(eta,u(3,:));
hold on;
end
end
% --------------------------------------------------------------------------
function dy=equation(~,y)
%dy=zeros(4,1);
dy=[y(2);
y(3);
-y(1)*y(3);
%---------------------------------------------------------------------
];
end
function res = bvpbc(y0,yinf,r_actual)
res = [y0(1)
y0(2)-1+ r_actual
yinf(2)- r_actual
];
end
function yini=guess(x,ir,sol)
global ix
ix = ix + 1;
if ir==1
yini(1)=sol(1);
yini(2)=sol(2);
yini(3)=sol(3);
else
yini(1)=sol(1,ix);
yini(2)=sol(2,ix);
yini(3)=sol(3,ix);
end
end
Anil Gautam
Anil Gautam el 7 de Nov. de 2018
thanku very much. thanks a lot

Iniciar sesión para comentar.

Categorías

Más información sobre Specifying Target for Graphics Output en Help Center y File Exchange.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by