You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
shooting method for coupled ODE
1 view (last 30 days)
Show older comments
dd
ggg
gg
13 Comments
J. Alex Lee
on 23 Apr 2020
In your code, you're setting E' to g' in line 83...but really you want the RHS of line 84 to be assigned to dy(7)...right?
Is your issue that you have a BC on g itself, but you aren't solving for g?
J. Alex Lee
on 24 Apr 2020
Please wrap your updated code in the code wrapper (use alt+enter).
What error are you seeing?
T K
on 24 Apr 2020
error
Undefined function or variable 'dg'.
Error in proj/projfun (line 69)
dy(2) = (E*f)+(E*x*df)+(f*df)+((dg)/(alfa*Pr*rho))-(df/x)-(f/(x*x));
Error in bvparguments (line 105)
testODE = ode(x1,y1,odeExtras{:});
Error in bvp4c (line 130)
bvparguments(solver_name,ode,bc,solinit,options,varargin);
Error in proj (line 37)
sol= bvp4c(@projfun,@projbc,solinit,options);
(general code is attached as m file )
function sol= proj
clc;clf;clear;
%y1=F, y2=F', y3=t, y4=dt, y5=C, y6=dC, y7=E,y8=E',y9=g
KB = 1.3807e-15;%Boltzman Constant
Mu = 1e-2;
K = 1e5;
Kp = 40e5;
dp = (100e-9)*1e2;
rho = 997.1/1000;
C = 4179e4;
alfa = K/(rho*C);
%rhonf = 3970/1000;
%cnf = 765e4;
B = KB/(3*pi*Mu*dp);
Bs = ((0.26*K)/(2*K+Kp))*(Mu/rho);
myLegend1 = {};
myLegend2 = {};
rr = [4 5 6 7];
%pp =[0.3 0.6 0.7 0.9];
%qq =[4.81 4.9 4.95 5];
for i = 1:numel(rr)
Pr = rr(i);
% gamma=pp(i);
%gammmafi=qq(i);
%Pr=7;
Tinf = 300;
Fiinf = 0.045;
gamma = 0.05546;
gammafi = 0.0077;
DT = gamma*Tinf;
Dfi = gammafi*Fiinf;
y0 = [1,0,1,0,1,0,1,0,1];
options =bvpset('stats','on','RelTol',1e-4);
%m = linspace(0.01,1);
m = linspace(0.1,0.44);
solinit = bvpinit(m,y0);
sol= bvp4c(@projfun,@projbc,solinit,options);
%y1=deval(sol,0)
solinit= sol;
figure(1)
plot(sol.x,(sol.y(9,:)))
grid on,hold on
myLegend1{i}=['pr = ',num2str(rr(i))];
figure(2)
plot(sol.x,sol.y(9,:))
grid on,hold on
myLegend2{i}=['Pr= ',num2str(rr(i))];
i=i+1;
end
figure(1)
legend(myLegend1)
hold on
figure(2)
legend(myLegend2)
function dy= projfun(x,y)
dy= zeros(9,1);
% alignComments
f = y(1);
df = y(2);
t = y(3);
dt = y(4);
c = y(5);
dc = y(6);
E = y(7);
dE = y(8);
g = y(9);
dy(1) = df;
dy(2) = (E*f)+(E*x*df)+(f*df)+((dg)/(alfa*Pr*rho))-(df/x)-(f/(x*x));
dy(3) = dt;
dy(4) = (Pr*dt*x*E)+(Pr*f*dt)-(dt/x);
dy(5) = dc;
dy(6) = (1/(B*Tinf))*((alfa*Pr*dc*f)+(alfa*Pr*E*x*dc)-(((DT*Bs*Fiinf)/(Tinf*Dfi))*((dt/x)+((Pr*dt*x*E)+(Pr*f*dt)-(dt/x)))))-(dc/x);
dy(7) = (-1/x)*(4*E+df+(f/x));
dy(8) = (4*E*E) +(E*x*((-1/x)*(4*E+df+(f/x)))+f*((-1/x)*(4*E+df+(f/x))) -(((-1/x)*(4*E+df+(f/x)))/x));
dy(9) = dg ;
end
end
function res= projbc(ya,yb)
res= [ya(1)-1;
ya(3)-1;
ya(5)-1;
ya(7)-1;
ya(9)-1;
yb(3);
yb(5);
yb(7);
yb(8)];
end
J. Alex Lee
on 25 Apr 2020
The error is straightforward...dg is undefined.
The reason is that you don't have equations in a form where all of the derivatives are explicitly written out (you don't have an explicit form dg/dx = ...)
It's not clear to me that it is possible to express your system in terms of 9 explicit derivatives.
If not, here's somet things to get started:
I don't think you would be able to use the "shooting" method either.
T K
on 26 Apr 2020
Edited: T K
on 26 Apr 2020
Welcome Doctor J.Alex Lee .
I am very happy and grateful for your effort with me.
I have expressed my system with 10 explicit derivatives As shown in line (77&78) in the code and there is no error appear when run the code.
From your point of view, is the code correct (by expressing the system with 10 explicit derivatives instead of 9 explicit derivatives) ????
T K
on 26 Apr 2020
function sol= proj
clc;clf;clear;
%y1=F, y2=F', y3=t, y4=dt, y5=C, y6=dC, y7=E,y8=E',y9=g
KB = 1.3807e-15;%Boltzman Constant
Mu = 1e-2;
K = 1e5;
Kp = 40e5;
dp = (100e-9)*1e2;
rho = 997.1/1000;
C = 4179e4;
alfa = K/(rho*C);
%rhonf = 3970/1000;
%cnf = 765e4;
B = KB/(3*pi*Mu*dp);
Bs = ((0.26*K)/(2*K+Kp))*(Mu/rho);
myLegend1 = {};
myLegend2 = {};
rr = [4 5 6 7];
%pp =[0.3 0.6 0.7 0.9];
%qq =[4.81 4.9 4.95 5];
for i = 1:numel(rr)
Pr = rr(i);
% gamma=pp(i);
%gammmafi=qq(i);
%Pr=7;
Tinf = 300;
Fiinf = 0.045;
gamma = 0.05546;
gammafi = 0.0077;
DT = gamma*Tinf;
Dfi = gammafi*Fiinf;
y0 = [1,0,1,0,1,0,1,0,1,0];
options =bvpset('stats','on','RelTol',1e-4);
%m = linspace(0.01,1);
m = linspace(0.008,0.2);
solinit = bvpinit(m,y0);
sol= bvp4c(@projfun,@projbc,solinit,options);
%y1=deval(sol,0)
solinit= sol;
figure(1)
plot(sol.x,(sol.y(9,:)))
grid on,hold on
myLegend1{i}=['pr = ',num2str(rr(i))];
figure(2)
plot(sol.x,sol.y(10,:))
grid on,hold on
myLegend2{i}=['Pr= ',num2str(rr(i))];
i=i+1;
end
figure(1)
legend(myLegend1)
hold on
figure(2)
legend(myLegend2)
function dy= projfun(x,y)
dy= zeros(9,1);
% alignComments
f = y(1);
df = y(2);
t = y(3);
dt = y(4);
c = y(5);
dc = y(6);
E = y(7);
dE = y(8);
g = y(9);
dg=y(10);
dy(1) = df;
dy(2) = (E*f)+(E*x*df)+(f*df)+((dg)/(alfa*Pr*rho))-(df/x)-(f/(x*x));
dy(3) = dt;
dy(4) = (Pr*dt*x*E)+(Pr*f*dt)-(dt/x);
dy(5) = dc;
dy(6) = (1/(B*Tinf))*((alfa*Pr*dc*f)+(alfa*Pr*E*x*dc)-(((DT*Bs*Fiinf)/(Tinf*Dfi))*((dt/x)+((Pr*dt*x*E)+(Pr*f*dt)-(dt/x)))))-(dc/x);
dy(7) = (-1/x)*(4*E+df+(f/x));
dy(8) = (4*E*E) +(E*x*((-1/x)*(4*E+df+(f/x)))+f*((-1/x)*(4*E+df+(f/x))) -(((-1/x)*(4*E+df+(f/x)))/x));
dy(9) = dg ;
dy(10)=0;
end
end
function res= projbc(ya,yb)
res= [ya(1)-1;
ya(3)-1;
ya(5)-1;
ya(7)-1;
ya(9)-1;
yb(1)-0.1;
yb(3);
yb(5);
yb(7);
yb(8)];
end
J. Alex Lee
on 26 Apr 2020
You are now just asserting that
, i.e., that g is constant (with value 1 according to your 5th BC). It's unlikely what you want. If it is, just remove g completely from the equations and insert 0 in place of dg in Eq. 2.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/287731/image.png)
I'm also unclear where the new BC is coming from (yb(1)-0.1=0) associated with the new DoF.
So no, I don't think this is what you want...
J. Alex Lee
on 27 Apr 2020
I do not understand what you mean.
J. Alex Lee
on 27 Apr 2020
You've asserted that g' does not vary in x by the line, i.e., g''=0
dy(10)=0;
which may or may not be true of your equations, I don't know.
So rather than g being constant (as I suggested before incorrectly), g' will be constant and take whatever value is consistent with your new boundary condition
yb(1)-0.1
Since I don't know the origin of this BC, I don't know if it is consistent with your original problem. I'm not used to dealing with implicit forms, maybe I'm missing something basic. Hope someone else can chime in at this point, I'm at my limit.
See Also
Categories
Find more on Graphics Performance in Help Center and File Exchange
Tags
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- 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)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)