Multi variable simultaneous differential equations with nested function

I would like to solve a set of mass transfer differential equations which are not only over time but over time and reactor length. I know I have to do a nested function but not quite sure how.

ds/dt=(d^2 s)/dx^2 - ds/dx- s

dq/dt=-∝.s

the 4 dimensionless variables are s (solvent concentration), x (reactor length) and t (time) and q (solute concentration). alpha is a constant.

The complexity is that while dq and ds are over time, there is a 2nd order d2s/dx2 term as well as ds/dx term in equation one, i.e. s is a factor of time AND distance inside reactor.

ANY help would be much much appreciated.

Answer 1

[John BG]

Hi Isaac

Modifying the initial example of the function pdepe please not your s is here u:

function pdepe_exercise_01

m = 0;
x = linspace(-10,10,200);
t = linspace(0,2,50);

sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);    % Extract first solution component as u
u = sol(:,:,1);

figure(1);surf(x,t,u)      % surface plot 
title('Numerical solution computed with 20 mesh points.')
xlabel('Distance x')
ylabel('Time t')

6
figure(2)   % A solution profile 
plot(x,u(end,:))
title('Solution at t = 2')
xlabel('Distance x')
ylabel('u(x,2)')

% --------------------------------------------------------------
function [c,f,s] = pdex1pde(x,t,u,DuDx)
c = 1;
f = DuDx;
s = -DuDx-u;

% --------------------------------------------------------------
function u0 = pdex1ic(x)
u0 = sin(pi*x);
% --------------------------------------------------------------

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = pi * exp(-t);
qr = 1;

setting c=pi^2

Since you neither defined initial or boundary conditions i left the same as in the example.

If you find this answer of any help solving your question,

please click on the thumbs-up vote link,

thanks in advance

John

[e-mail]

[Isaac Bt]

Hi John,

Really appreciate your prompt response! Thrilled!

You're right, I missed the initial conditions cause there are a lot of constants which I'd removed for simplicity. Some initial & boundary conditions are

t=0,s=0;

x=0,s=1;

x=1,ds/dx=0;

t=0,q=0;

So you don't think a nested function is required here?

What I need to get from these governing equations is plots like below where (pi= s) vs (tau =t) and q vs. x (solute conc.'n vs dimensionless distance).

Thanks again for being so great and for your time.

[John BG]

pdepe initial conditions go in the function @pdex1ic modify this function accordingly.

pdepe boundary conditions go in function @pdex1bc

Check the MATLAB help for pdepe, you will see more useful examples.

[Isaac Bt]

John, Aren't my governing equations ODEs? Looking at the pdepe example I guess I should be able to make it fit my equations.

[John BG]

from MATLAB help PDEPE stands for partial differential Equations with parabolic / elliptical solutions. The kind that solutions have to do with

exp(j* ..)

for instance Electromagnetic fields trapped in a cavity, you have the waves varying with time, and [x y z]

ODE stands for ordinary differential equation: the equations have the same

f(DN_y  Dn-1_y , .. , y''', y'', y', y, g(t))

basic intro: http://www.math.vanderbilt.edu/~daphne/courses/math226/MatlabODEs.pdf

ODE systems would be

f1(DN_y1  Dn-1_y1 , .. , y1''', y1'', y1', y1, g1(t))
f(DN_y2  Dn-1_y2 , .. , y2''', y2'', y2', y2, g2(t))
f(DN_y3  Dn-1_y3 , .. , y3''', y3'', y3', y3, g3(t))

as complex as ODE system may be their initial conditions and boundary conditions are relatively simple compare to those needed in PDE.

Earlier on I did not mention your second equation because that is and ODE: you solve the first one pdepe, and the you simply get q time integrating s.

If you need to build an ODE PDE knowledge base mathworks website, MATLAB help and MIT Openware are great, among many possible references

Regards

John

[e-mail]

[RB]

Hello! I am trying to solve the following system of ODEs I have the terminal values but the equation is involving 2*2 matrices. Solve the following system of ODE's

d/dt Y(t)=-(H(t))^{'}Y(t) d/dt V(t)=-Y(t)^{'}Q^{'}QY(t)

with the terminal conditions Y(T)=I_{2} and V(T)=0. Y, V, H and Q are 2*2 matrices and ' stands for transpose.

Also H(t)=H-Q^{'}Q*W(T-t) where W(T-t)=W(tau) is determined by the matrix Ricatti equation

dW/dtau = W(tau)H+H^{'}W(tau)-2W(tau)Q^{'}QW(tau)+I_{2} such that at tau=0 W(0)=0. I am badly stuck up. I am unable to use ode45 since I want matrix result. It would be great if you can give a suggestion. \[ H= \begin{bmatrix} -0.5 & 0.4 \\ 0.007 & -0.008 \end{bmatrix} \] \[ Q= \begin{bmatrix} 0.06 & -0.0006 \\ -0.06 & 0.006 \end{bmatrix} \] Thanks. <http://quicklatex.com/cache3/ef/ql_28653ec4bfcd28208120ea76b9009bef_l3.png>I have tried everything in MATLAB but I can not get a solution. Any kind of help would be great.

Actually I am solving a system of matrix odes where dx/dt, dy/dt, dz/dt are all 2 by 2 matrices and I need to get X,Y,Z which are again 2 by 2 matrices. My code is function [dXdt, dYdt, dVdt] = mwish7(t, X, Y, V, A, B, R) X = reshape(X, size(A)); Y = reshape(Y, size(A)); V = reshape(V, size(A)); dXdt = -A.'*X - X*A + 2*X*(B.'*B)*X - R; dYdt = -(A-(B.'*B)*X).'*Y; dVdt = -Y.'*(B.'*B)*Y; dXdt = dXdt(:); dYdt = dYdt(:); dVdt = dVdt(:); end and the ode45 program I use is A=[-0.5, 0.4; 0.007, -0.008]; B=[0.06 -0.0006; -0.06, 0.006]; R = [1 0; 0 1]; X0 = [0, 0; 0, 0]; Y0 = [1 0; 0 1]; V0 = [0, 0; 0, 0]; [T X Y V] = ode45(@(t,X,Y,V)mwish7(t, X, Y, V, A, B, R), [15 0], X0, Y0, V0) I get the error: NOT ENOUGH INPUT ARGUMENTS.

Thanks.

Regards, RB