I am not getting multiple graphs(iterative) when I run the code for a coupled bvp ODE using bvp4c
2 次查看(过去 30 天)
显示 更早的评论
naygarp
2017-12-5
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
采纳的回答
Torsten
2017-12-6
https://de.mathworks.com/help/matlab/ref/hold.html
Best wishes
Torsten.
21 个评论
naygarp
2017-12-6
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
2017-12-7
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
2017-12-9
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
2017-12-11
If the code you posted worked, reset the parameters alp, lambda and pp (resp. gama) to its earlier values.
Best wishes
Torsten.
naygarp
2017-12-11
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
2017-12-12
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/187365/image.jpeg)
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
2017-12-12
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
2017-12-12
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
2017-12-12
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
2017-12-12
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
2017-12-13
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
2017-12-13
Yes, it's really alarming how much erroneous results are published under the pressure to "publish or perish".
Best wishes
Torsten.
naygarp
2017-12-14
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.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Visual Exploration 的更多信息
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 (한국어)