이 질문을 팔로우합니다.
- 팔로우하는 게시물 피드에서 업데이트를 확인할 수 있습니다.
- 정보 수신 기본 설정에 따라 이메일을 받을 수 있습니다.
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.
naygarp
2017년 12월 13일
The nature of your graphs are similar to what I got with the help of your code
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 (한국어)