Solving a coupled system of ODE's that depend on different variables
1 view (last 30 days)
Show older comments
I have a set of coupled ODE's that I need to solve. However, the ODE's are not dependent on the same variable. For example I have the following first order system:
dA0/dx = A,
dA/dx = f(A0(x), A(x), F(y), y)
dF/dy = f(A0(x), A(x), y).
For example, the dF/dy looks like: p*(y-y^2)*(A0*A+A^2) where p is a constant And the dA/dx equation looks like: ([p*A0*(y-y^2)-A0]*A)/(F+A*(y-y^2))
The first dA0/dx is to get them all down to first order.
Is there a numerical or analytic method to solve these since the derivatives are dependent on different variables? All the solutions I can find need dA and dF to be with respect to the same variable and I can't get ODE45 to work with these equations.
Here is my current code:
dx=(L/H)/100; x=(0:dx:(L/H))'; y=[0:H/(size(x,1)-1):H]; %Pressure drop through Unexposed area %Pressure drop count=1; % Equation: A3 = (420/(11*Re0))+(g(i)*h(i)+84/11*f(i)*g(i)-35/11*NClam/Re0*h(i)-70/11*Cturb*g(i)*h(i)-70/11*1/Re0*h(i))/f(i); while count <2 f=zeros(length(x),1); % f=A g=zeros(length(x),1); % g=A' h=zeros(length(x),1); % h=A'' A3=zeros(length(x),1); % A3=A''' f(1)=-.99; %Boundary Condition at dead end g(1)=0; %Boundary Condition at dead end %1st 2 Shooting Guesses SG(1)=.0001; SG(2)=.00012;
for j=1:2
h(1)=SG(j);
for i=1:size(x,1)-1
% Runge Kutta Slopes
k1f=dx*g(i);
k1g=dx*h(i);
k1h=dx*((420/(11*Re0))+(g(i)*h(i)+84/11*(1+f(i))*g(i)-35/11*NClam/Re0*h(i)-70/11*Cturb*g(i)*h(i)-70/11*1/Re0*h(i))/(1+f(i)));
k2f=dx*(g(i)+1/2*k1g);
k2g=dx*(h(i)+1/2*k1h);
k2h=dx*((420/(11*Re0))+((g(i)+1/2*k1g)*(h(i)+1/2*k1h)+84/11*(1+(f(i)+1/2*k1f))*(g(i)+1/2*k1g)-35/11*NClam/Re0*(h(i)+1/2*k1h)-70/11*Cturb*(g(i)+1/2*k1g)*(h(i)+1/2*k1h)-70/11*1/Re0*(h(i)+1/2*k1h))/(1+(f(i)+1/2*k1f)));
k3f=dx*(g(i)+1/2*k2g);
k3g=dx*(h(i)+1/2*k2h);
k3h=dx*((420/(11*Re0))+((g(i)+1/2*k2g)*(h(i)+1/2*k2h)+84/11*(1+(f(i)+1/2*k2f))*(g(i)+1/2*k2g)-35/11*NClam/Re0*(h(i)+1/2*k2h)-70/11*Cturb*(g(i)+1/2*k2g)*(h(i)+1/2*k2h)-70/11*1/Re0*(h(i)+1/2*k2h))/(1+(f(i)+1/2*k2f)));
k4f=dx*(g(i)+k3g);
k4g=dx*(h(i)+k3h);
k4h=dx*((420/(11*Re0))+((g(i)+k3g)*(h(i)+k3h)+84/11*(1+(f(i)+k3f))*(g(i)+k3g)-35/11*NClam/Re0*(h(i)+k3h)-70/11*Cturb*(g(i)+k3g)*(h(i)+k3h)-70/11*1/Re0*(h(i)+k3h))/(1+(f(i)+k3f)));
% Evaluation
f(i+1)=f(i)+1/6*(k1f+2*k2f+2*k3f+k4f);
g(i+1)=g(i)+1/6*(k1g+2*k2g+2*k3g+k4g);
h(i+1)=h(i)+1/6*(k1h+2*k2h+2*k3h+k4h);
end
% R is f at channel end, we want equal to Uout
R(j)=f(size(x,1));
err=abs(uout-(R(j)));
end
% This loop iterates to find the right guess for f''(0) so that
% f(x/L)= Uout
while err>1e-9
j=j+1;
% New guess is found by linear interposation of the last 2 guesses
SG(j)=SG(j-2)+(SG(j-1)-SG(j-2))/(R(j-1)-R(j-2))*(-R(j-2));
h(1)=SG(j);
for i=1:size(x,1)-1
% Runge Kutta Slopes
k1f=dx*g(i);
k1g=dx*h(i);
k1h=dx*((420/(11*Re0))+(g(i)*h(i)+84/11*(1+f(i))*g(i)-35/11*NClam/Re0*h(i)-70/11*Cturb*g(i)*h(i)-70/11*1/Re0*h(i))/(1+f(i)));
k2f=dx*(g(i)+1/2*k1g);
k2g=dx*(h(i)+1/2*k1h);
k2h=dx*((420/(11*Re0))+((g(i)+1/2*k1g)*(h(i)+1/2*k1h)+84/11*(1+(f(i)+1/2*k1f))*(g(i)+1/2*k1g)-35/11*NClam/Re0*(h(i)+1/2*k1h)-70/11*Cturb*(g(i)+1/2*k1g)*(h(i)+1/2*k1h)-70/11*1/Re0*(h(i)+1/2*k1h))/(1+(f(i)+1/2*k1f)));
k3f=dx*(g(i)+1/2*k2g);
k3g=dx*(h(i)+1/2*k2h);
k3h=dx*((420/(11*Re0))+((g(i)+1/2*k2g)*(h(i)+1/2*k2h)+84/11*(1+(f(i)+1/2*k2f))*(g(i)+1/2*k2g)-35/11*NClam/Re0*(h(i)+1/2*k2h)-70/11*Cturb*(g(i)+1/2*k2g)*(h(i)+1/2*k2h)-70/11*1/Re0*(h(i)+1/2*k2h))/(1+(f(i)+1/2*k2f)));
k4f=dx*(g(i)+k3g);
k4g=dx*(h(i)+k3h);
k4h=dx*((420/(11*Re0))+((g(i)+k3g)*(h(i)+k3h)+84/11*(1+(f(i)+k3f))*(g(i)+k3g)-35/11*NClam/Re0*(h(i)+k3h)-70/11*Cturb*(g(i)+k3g)*(h(i)+k3h)-70/11*1/Re0*(h(i)+k3h))/(1+(f(i)+k3f)));
% Evaluation
f(i+1)=f(i)+1/6*(k1f+2*k2f+2*k3f+k4f);
g(i+1)=g(i)+1/6*(k1g+2*k2g+2*k3g+k4g);
h(i+1)=h(i)+1/6*(k1h+2*k2h+2*k3h+k4h);
end
% R is f at channel end, we want equal to Uout
R(j)=f(size(x,1));
err=abs(uout-(R(j)));
end
% Calculate A'''
for i = 1:length(x)
A3(i) = (420/(11*Re0))+(g(i)*h(i)+84/11*(1+f(i))*g(i)-35/11*NClam/Re0*h(i)-70/11*Cturb*g(i)*h(i)-70/11*1/Re0*h(i))/(1+f(i));
end
% Calculate U and V
for n = 1:length(x)
for m = 1:length(y)
U(n,m) = uout.*(1+f(n)).*(6*(y(m)/H).*(1-(y(m)/H)));
V(n,m) = -uout.*g(n)*(3*(y(m)/H).^2-2*(y(m)/H).^3);
end
end
count = count+1;
end
3 Comments
ragesh r menon
on 26 Jun 2014
I use ode45 to solve coupled nonlinear ODEs and get accurate results. For example my state equations are
function dxdt=LOS(t,x)
......
dxdt(1)=Vs*cos(x(5)-x(2))-Vt*cos(x(7)-x(2));%LOS separation
dxdt(2)=(Vs*sin(x(5)-x(2))-Vt*sin(x(7)-x(2)))/x(1);%sLOS rate
......
end
Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!