How do I use MatLab to solve this set of differential equations?
3 views (last 30 days)
Show older comments
CarlosRobelli
on 27 Mar 2015
Commented: Star Strider
on 28 Mar 2015
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?
0 Comments
Accepted Answer
Star Strider
on 27 Mar 2015
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.
4 Comments
Star Strider
on 28 Mar 2015
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.)
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!