How do I use MatLab to solve this set of differential equations?

So here is my code (The code only contains the initial constants and equations):

%Constants
Fb = .103;
k = 40/60*1000;
alpha = .010235;
mass = 30;
vo = 8.571;
epsilon = .6;

%Equations defining the Concentrations and the Rate of Destruction
Ca = @(x,y) Fb*(4-2*x)*y/((1+epsilon*x)*vo);
Cb = @(x,y) Fb*(1-x)*y/((1+epsilon*x)*vo);
rb = @(x,y) k*Ca(x,y)^2*Cb(x,y)^2;

% Differential Equations
dx = @(x,y) rb(x,y)/Fb;
dy = @(x,y) -alpha*(1+epsilon*x)/(2*y);

x is the percent conversion and it is a function of w (catalyst mass). y is a ratio of pressures and is also a function of w. When w = 0, x = 0 and y = 1. I need to find the percent conversion (x) for a given catalyst mass (w).

dx represents dx/dw and dy represents dy/dw.

How do I solve a problem like this using MatLab?

Answer 1

[Star Strider]

I’m not certain what your initial conditions are, so I used x=0, y=1. Leaving the rest of your code unchanged (constants and function defined before the code here), the differential equation integration is:

%Constants
Fb = .103;
k = 40/60*1000;
alpha = .010235;
mass = 30;
vo = 8.571;
epsilon = .6;

Ca = @(x,y) Fb*(4-2*x)*y/((1+epsilon*x)*vo);
Cb = @(x,y) Fb*(1-x)*y/((1+epsilon*x)*vo);
rb = @(x,y) k*Ca(x,y)^2*Cb(x,y)^2;

% MAPPING: xy(1) = x,  xy(2) = y
DEs = @(t,xy) [rb(xy(1),xy(2))/Fb;  -alpha*(1+epsilon*xy(1))./(2*xy(2))];

tspan = linspace(0, 94.9, 100);
sy0 = [0; 1];
[t,xy] = ode15s(DEs, tspan, sy0);

figure(1)
plot(t, xy)
grid
legend('x(t)', 'y(t)')

There is a singularity or some sort of discontinuity at about t=94.9, so I only took ‘tspan’ out to that limit. I know nothing about your system, so I just left it at that. I don’t see ‘w’ in your equations anywhere, so I will leave that to you. I am simply showing you how to integrate them in MATLAB.

[CarlosRobelli]

Thanks! Sorry about the confusion of the 'w' variable. In my code dx and dy are actually meant to represent dx/dw and dy/dw. Where both x and y are functions of w.

I hope that clarifies things.

Would that change your solution?

[Star Strider]

My pleasure!

It definitely clarifies things, and it means that the initial conditions I defined are correct.

The only change in my code would be to change ‘t’ to ‘w’, so it becomes:

% MAPPING: xy(1) = x,  xy(2) = y
DEs = @(w,xy) [rb(xy(1),xy(2))/Fb;  -alpha*(1+epsilon*xy(1))./(2*xy(2))];

wspan = linspace(0, 94.9, 100);
xy0 = [0; 1];
[w,xy] = ode15s(DEs, wspan, xy0);

figure(1)
plot(w, xy)
grid
legend('x(w)', 'y(w)')

Does it do what you want?

[CarlosRobelli]

It does, thanks!

I have another question if you don't mind. Let's say that I want to set w = 20 and I want to see what the values for x and y would be at that point. How would I find these values? Would I use the equations defined in this line:

DEs = @(w,xy) [rb(xy(1),xy(2))/Fb;  -alpha*(1+epsilon*xy(1))./(2*xy(2))];

[Star Strider]

My pleasure!

That is the correct line, since it defines your differential equation system. You won’t need it, though, other than to use it to integrate your system.

Setting ‘w’ to 20 and finding the ‘x’ and ‘y’ values there is easy. Do the integration as I outlined. Then use the interp1 function on ‘w’ and ‘xy’ that are the solutions to your equation, and interpolate:

xy20 = interp1(w, xy, 20);

The ‘xy20’ variable will have the values of ‘x’ and ‘y’ at w=20. Check them with the plot.

(Note: I didn’t actually run the interp1 code with your integrated system, but I’ve enough experience with interp1 to be confident that it will work. If you encounter any problems with it, let me know.)