I am not getting multiple graphs(iterative) when I run the code for a coupled bvp ODE using bvp4c

Two ODEs are:

F"=G(G^2+gamma^2)/(G^2+lambda*gamma^2)

G'= 1/3F'^2-2/3(F*F")

subject to: F(xi)=alpha/2, F'(xi)=1 at xi=0 where 'alpha' is a parameter (wall parameter)

F'(xi)= 0 as xi tends to infinity

I should be getting a multiple graphs varying the parameter ' alpha'

The code that I have run is:

function sol= proj
clc;clf;clear;
global lambda gama alp
lambda=0.5;
gama=1;
pp=[0:0.5:1.0];
figure(1)
plot(2,1);
solinit= bvpinit([0:0.01:10],[0,1,0]);
for alp= pp
sol= bvp4c(@projfun,@projbc,solinit);
solinit= sol;
plot(2,1);plot(sol.x,sol.y(2,:))
end
end
function f= projfun(x,y)
global lambda gama
f= [y(2);y(3)*(y(3)^2+gama.^2)/(y(3)^2+lambda*gama.^2);y(2)^2/3-(2*y(1)*y(3)*(y(3)^2+gama.^2))/(3*(y(3)^2+lambda*gama.^2))];
end
function res= projbc(ya,yb)
global alp
res= [ya(1)-alp/2; ya(2)-1.0; yb(2)];
end

Answer 1

[Torsten]

https://de.mathworks.com/help/matlab/ref/hold.html

Best wishes

Torsten.

[naygarp]

Thanks once again.

[naygarp]

I am attaching a research paper from where the above equations are taken. I could not figure out how the graphs illustrated in fig. 4 (a,b,c). where f''(0) are plotted against gamma are obtained. Could you please help.

[Torsten]

pp=[0 1 2 3 4];
lambda=0.5;
alp=0.5;
for i=1:numel(pp)
  gama=pp(i);
  sol= bvp4c(@projfun,@projbc,solinit);
  y=deval(sol,0);
  G=y(3);
  Fss(i)=G*(G^2+gama^2)/(G^2+lambda*gama^2);
  solinit=sol;
end
plot(pp,Fss)

Best wishes

Torsten.

[naygarp]

I have integrated the above code into my earlier code and tried to run but, the following error occurs:

Error using bvp4c (line 251) Unable to solve the collocation equations -- a singular Jacobian encountered.

Error in proj_variable_thickness (line 12) sol= bvp4c(@projfun,@projbc,solinit);

 function sol= proj
clc;clf;clear;
global lambda gama alp
alp=0.5;
lambda=0.5;
pp=[0 1 2 3 4];
figure(1)
solinit= bvpinit([0:0.01:10],[0,1,0]);
for i=1:numel(pp)
    gama=pp(i);
sol= bvp4c(@projfun,@projbc,solinit);
solinit=sol;
y=deval(sol,0);
G=y(3);
Fss(i)=G*(G^2+gama^2)/(G^2+lambda*gama^2);
plot(pp,Fss)
end
end
function f= projfun(x,y)
global lambda gama
f= [y(2);y(3)*(y(3)^2+gama^2)/(y(3)^2+lambda*gama^2);y(2)^2/3-(2*y(1)*y(3)*(y(3)^2+gama^2))/(3*(y(3)^2+lambda*gama^2))];
end
function res= projbc(ya,yb)
global alp
res= [ya(1)-alp/2; ya(2)-1.0; yb(2)];
end

[Torsten]

If the code you posted worked, reset the parameters alp, lambda and pp (resp. gama) to its earlier values.

Best wishes

Torsten.

[naygarp]

No, the code didn't run, there's two errors

 Error using bvp4c (line 251) Unable to solve the collocation equations -- a singular Jacobian encountered.
Error in proj_variable_thickness (line 12) sol= bvp4c(@projfun,@projbc,solinit);

[naygarp]

I have used 'deval' and got the values of y(3) at xi=0 and then used the value in the equation as you suggested

F" = y(3)*(y(3)^2+gama^2)/(y(3)^2+lambda*gama^2)

then I used the following code to get the graphs But the graphs obtained do not match exactly but the nature of the graphs are similar as given in the research paper I posted earlier.

Please help where have I gone wrong.

202.jpg>>

>>

I used a different set of codes to get those graphs

clc; clf; clear;
global lambda gama
lambda=0.5;
qq=[-0.5111 -0.6150 -0.8568]
pp=[0.2:0.2:8];
for i=1:numel(pp);
  k = 1;
  for G=qq;    
      gama=pp(i); 
      Fss(i,k)= G*(G^2+gama^2)/(G^2+lambda*gama^2);
      k = k + 1;
  end
end
plot(pp,Fss);hold on

[Torsten]

You will get different values for G for different values of gama.

So the length of "qq" must be the same as the length of "pp" in your code.

Assuming that you kept "alp" and "lambda" constant, you will only get 1 curve as graph, not 3.

Best wishes

Torsten.

[naygarp]

The values of G for 3 different values of 'alpha' are obtained I guess as shown in the research paper. That's I think how those 3 curves are obtained.

[Torsten]

Let's take figure 4a).

The curve for alpha=0.5 (the green curve), e.g., is obtained as follows:

Fix alpha=0.5 and lambda=0.5 in the ODE model.

Vary gamma in the range 0-4, lets say gamma=[0 1 2 3 4].

Then, for each gamma from the list, solve the ODE using bvp4c.

Use "deval" to evaluate F''(0) for all values of gamma.

This will also give 5 values for F''(0).

Then plot these values for F''(0) against the gamma vector.

This will give the green curve in figure 4a).

But I already gave you the code for this plot. Why don't you use it ?

Best wishes

Torsten.

[naygarp]

Thanks again for making it clear.

I have used the code as you stated in the ODE model, but I am receiving some error, it's not running Please check, where I have done the mistake:

function sol= proj
clc;clf;clear;
global lambda gama alp
pp=[0 1 2 3 4];
alp=0.5;
lambda=0.5;
figure(1)
solinit= bvpinit([0:0.01:5],[0,1,0]);
for i=1:numel(pp)
  gama=pp(i)
sol= bvp4c(@projfun,@projbc,solinit);
y=deval(sol,0)
G= y(3);
Fss(i)=G*(G^2+gama^2)/(G^2+lambda*gama^2);
solinit=sol;
end
plot(qq,Fss)
end
function f= projfun(x,y)
global lambda gama
f= [y(2);y(3)*(y(3)^2+gama^2)/(y(3)^2+lambda*gama^2);y(2)^2/3-(2*y(1)*y(3)*(y(3)^2+gama^2))/(3*(y(3)^2+lambda*gama^2))];
end
function res= projbc(ya,yb)
  global alp
res= [ya(1)-alp/2; ya(2)-1.0; yb(2)];
end

[naygarp]

The code runs for

pp=[0.2 1 2 3 4 5]

But the graph I obtained is dissimilar,

[Torsten]

And what curve do you get for lambda = 2 ?

Best wishes

Torsten.

[naygarp]

This is for lambda =2

[naygarp]

This is for lambda =1

[Torsten]

So it's almost constant for lambda = 1.

Now if the curve was decreasing for lambda = 0.5, your curves look quite similar to the curves in the article, don't they ?

Best wishes

Torsten.

[naygarp]

In the research paper the graph shown is actually increasing for lambda =0.5 isn't it and decreasing for lambda= 2.

What I obtained is opposite or maybe the details in the paper could be wrong.

[Torsten]

That's what I get with COMSOL:

[naygarp]

The nature of your graphs are similar to what I got with the help of your code

[Torsten]

Yes, it's really alarming how much erroneous results are published under the pressure to "publish or perish".

Best wishes

Torsten.

[naygarp]

Yes after observing two research papers from well-known publishers, it's heartening to see how such details get missed. Anyway, thanks a lot for guiding me along.