Boundary condition in non-linear ODE

Dear all,

I wanted to solve this two set of non-linear ODE using matlab :

are constant

The boundary conditions are the following :

and at

and at ( ν is an arbitrary constant < 1)

this the code that I constracted so far

function bvp4c_mathworks
rspan = [0.01 1];
init = zeros(1,4);
solinit = bvpinit(rspan,init);
sol = bvp4c(@ode4,@bc4,solinit);
eta = sol.x;
theta = sol.y(1,:);
Sr  = sol.y(2,:);
plot(eta,theta)
hold on
plot(eta,Sr,'r')
hold off
legend('Nr(r)','\beta(r)')
end
function du = ode4(eta,u)
theta  = u(1);
Sr   = u(2);             % beta
dtheta = u(3);             % d(theta)/dr
dSr  = u(4);             % d(Sr)/dr
lambda =15.94;
P=12;                    %P=F*a/D; F is the applied force ; a radius of the membrane ; D = E*h^3/12(1-nu^2)
alpha = 3;               %alpha =C*a^2/D ; C in-plane stifnnes
du(1) = dtheta;
du(2) = dSr;
du(3) = (P/(2*pi*eta)-(1/eta)*dtheta+(1/eta^2+lambda^2+Sr));
du(4) = (alpha*theta^2/(2*eta^2)+3/eta*dSr);
du(4) = du(4)/eta;
end
function res = bc4(u0, ur)
res = [ur(1)-0
       ur(2)-0
       ur(3)-0
       u0(2)-0];
end

[Torsten]

du(3) and du(4) and your boundary conditions do not correspond to your mathematical equations.

[zakaria azdad]

I put theta= u(3) than du(3) = d²theta, probably I am wrong, could you give a suggestion to fix the bugg? and how could you modify the boundary conditions.

thanks

Answer 1

[Torsten]

function du = ode4(eta,u)
  theta  = u(1);
  Sr     = u(2);             % beta
  dtheta = u(3);             % d(theta)/dr
  dSr    = u(4);             % d(Sr)/dr
  lambda = 15.94;
  P=12;                      %P=F*a/D; F is the applied force ; a radius of the membrane ; D = E*h^3/12(1-nu^2)
  alpha = 3;                 %alpha =C*a^2/D ; C in-plane stifnnes
  du = zeros(4,1);
  du(1) = dtheta;
  du(2) = dSr;
  du(3) = P/(2*pi*eta)- 1/eta*dtheta + theta*(1/eta^2+lambda^2+Sr);
  du(4) = -alpha*theta^2/(2*eta^2) - 3/eta*dSr;
end

function res = bc4(ul,ur)
  nu = 0.1;
  res = zeros(4,1);
  res(1) = ul(1);
  res(2) = ur(1);
  res(3) = ul(4);
  res(4) = ur(4)+(1-nu)*ur(2);
end

[zakaria azdad]

I wanted to solve for different value of P and store the output of theta and Sr for each value of P.

function du = ode4(eta,u)
  theta  = u(1);
  Sr     = u(2);             % beta
  dtheta = u(3);             % d(theta)/dr
  dSr    = u(4);             % d(Sr)/dr
  lambda = 15.94;
  P=linspace(1e-3,1e3,1e4);  %P=F*a/D; F is the applied force ; a radius of the membrane ; D = E*h^3/12(1-nu^2)
  alpha = 3;                 %alpha =C*a^2/D ; C in-plane stifnnes
  for i=1:length(P)
      du = zeros(4,i);
      du(1,i) = dtheta;
      du(2,i) = dSr;
      du(3,i) = P(i)/(2*pi*eta)- 1/eta*dtheta + theta*(1/eta^2+lambda^2+Sr);
      du(4,i) = -alpha*theta^2/(2*eta^2) - 3/eta*dSr;
  end
end

Than define so is there a guide how to do it?

[Torsten]

The loop over P must be in the main program "bvp4c_mathworks", not in ode4.

To obtain w, introduce a third differential equation

dw/dr = theta 

with initial condition

w(0) = 0

function bvp4c_mathworks
  rspan = [0.01 1];
  init = zeros(1,5);
  solinit = bvpinit(rspan,init);
  P = linspace(1e-3,1e3,1e4);  
  for i=1:numel(P)
    sol{i} = bvp4c(@(eta,u)ode4(eta,u,P(i),@bc4,solinit);
  end
end

function du = ode4(eta,u,P)
  theta  = u(1);
  Sr     = u(2);             % beta
  dtheta = u(3);             % d(theta)/dr
  dSr    = u(4);             % d(Sr)/dr
  lambda = 15.94;
  alpha = 3;                 %alpha =C*a^2/D ; C in-plane stifnnes
  du = zeros(5,1);
  du(1) = dtheta;
  du(2) = dSr;
  du(3) = P/(2*pi*eta)- 1/eta*dtheta + theta*(1/eta^2+lambda^2+Sr);
  du(4) = -alpha*theta^2/(2*eta^2) - 3/eta*dSr;
  du(5) = theta;
end  

function res = bc4(ul,ur)
  nu = 0.1;
  res = zeros(5,1);
  res(1) = ul(1);
  res(2) = ur(1);
  res(3) = ul(4);
  res(4) = ur(4)+(1-nu)*ur(2);
  res(5) = ul(5)
end

Best wishes

Torsten.

[zakaria azdad]

Thank you for your suggestions!!

I have modified the function so I can be able to plot different quantities in this way.

  rspan = linspace(0.01,1,100); 
  init = zeros(1,5);
  solinit = bvpinit(rspan,init);
  P = linspace(10,50,100);  
  for i=1:numel(P)
    sol{i} = bvp4c(@(eta,u)ode4(eta,u,P(i)),@bc4,solinit);
    eta{i} = sol{i}.x;
    theta{i} = sol{i}.y(1,:);
    Sr{i}  = sol{i}.y(2,:);
    W{i}   =sol{i}.y(5,:);
    plot(eta{1,i},P(i))
    hold on
  end

I tried to plot as a functionof P(i) but it doesn't work. do you see a an issue there?

I noticed that eta, theta and Sr have a structure of 1x100 cell and each cell have 45 colum. When I try to plot(w,P) it tells me that Vectors must be the same length.

However, when I change P = linspace(10,50,45); it works. How can I solve this issue?

[Torsten]

  rspan = linspace(0.01,1,100); 
  init = zeros(1,5);
  solinit = bvpinit(rspan,init);

  P = linspace(10,50,100);
  
  theta = zeros(numel(P),numel(rspan));
  Sr = zeros(numel(P),numel(rspan));
  w = zeros(numel(P));

  for i=1:numel(P)
    sol = bvp4c(@(eta,u)ode4(eta,u,P(i)),@bc4,solinit);
    theta(i,:) = deval(sol,rspan,1);
    Sr(i,:)  = deval(sol,rspan,2);
    w(i) = deval(sol,rspan(end),5);
  end

  plot(P,w)

[zakaria azdad]

I am bit lost.

Using your last suggestion, i can't determine eta vs theta. When I put eta =dval(sol,rspan) and than I plot etaVs theta, I get something different from using eta =sol.x and theta= sol.y(1,:) ?

can you clarify this issue?

thanks in advance.

[Torsten]

The theta for which P-value do you want to plot against eta ? For all of them ?

Then just use

plot(rspan,theta)

after the for-loop.

[zakaria azdad]

Hi, thaks a lot!!

Finaly I have the code put together. However, I get a warning message "Unable to meet the tolerance without using more than 2000 mesh points." eventhough my rspan is

rspan = linspace(0.001,1,5000);

I played arround with vlues and it just didn't converged.

This happened when a large value (E+7) of α are entred.

Any suggestions?

[Torsten]

Yes, use the solution of a converging run as initial guess for a subsequent run.

But you will have to do this on your own now because it's time to start learning MATLAB.

[zakaria azdad]

Any hint where to start?