Matlab ODE Solution

1
Matlab ODE Solution

• Matlab codes
• Matlab tutorial
• CSTR exercise

2
Matlab ODE Solution Codes
3
Matlab Tutorial

• Solution of ODE systems with Matlab
• ode23, ode45, ode113, ode15s, ode23s, ode23t,
  ode23tb
• Matlab functions for solving initial value
  problems for ordinary differential equations
• Syntax x,y solver(odefun,xspan,y0,options)
• solver one of the Matlab ODE solvers
• odefun name of function that evaluates the RHS
  of dy/dxf(x,y)
• xspan vector specifying the integration
  interval x0, xf
• y0 vector of initial conditions
• options solver options (optional)
• ode15i
• Matlab function for solving fully implicit
  differential equations
• See help ode15i for details

4
Matlab Tutorial cont.

• van der Pol equation as an ODE system
• Built-in Matlab functions for the van der Pol
  ODEs
• m1 vdp1(t,y)
• m1000 vdp1000(t,y)
• Use Matlab ODE solvers to find the solution
• Value of m effects the stiffness of the system
• small m, nonstiff, ode45 more efficient
• large m, stiff, stiff solver required

gtgt type vdp1
gtgt t,yode45(_at_vdp1,0 20,2 0)
plot(t,y(,1),'-o') gtgt t,yode15s(_at_vdp1000,0
3000,2 0) plot(t,y(,1),'-o')
gtgt tic, for i110 t,yode45(_at_vdp1,0 20,2
0) end toc gtgt tic, for i110
t,yode15s(_at_vdp1,0 20,2 0) end toc
Too slow, abort with Ctrl-C
gtgt t,yode45(_at_vdp1000,0 3000,2 0)
5
Nonisothermal Chemical Reactor

• Reaction A ? B
• Assumptions
• Pure A in feed
• Perfect mixing
• Negligible heat losses
• Constant properties (r, Cp, DH, U)
• Constant cooling jacket temperature

6
Matlab Exercise

• Parameter values
• k0 3.493x107 h-1, E 11843 kcal/kmol
• (-DH) 5960 kcal/kmol, rCp 500 kcal/m3/K,
• UA 150 kcal/h/K, R 1.987 kcal/kmol/K
• V 1 m3, q 1 m3/h,
• CAi 10 kmol/m3, Ti 298 K, Tc 298 K
• Problem
• Find the three steady-state points
• Integrate the model equations from each steady
  state
• Select one stable steady-state point
• Integrate the model equations from this initial
  condition for a different value of the coolant
  temperature Tc
• Plot the temperature response T(t)

7
Matlab Exercise Solution

• Find the three steady-state points
• Integrate the model equations from each steady
  state
• Select one stable steady-state point and
  integrate the model equations from this initial
  condition for a different value of the coolant
  temperature Tc
• Plot the temperature response T(t)

gtgt Tc 298 gtgt y1 fsolve(_at_cstr,9 300,,Tc)
y2fsolve(_at_cstr,5 330,,Tc) gtgt y3
fsolve(_at_cstr,2 360,,Tc) xo y1 y2
y3 xo 8.5637 311.1702 5.5176
339.1001 2.3591 368.0615
gtgt xo round(xo) gtgt df _at_(t,x,Tc)
cstr(x,Tc) gtgt ode45(df,0 100,xo(1,),,Tc) gtgt
ode45(df,0 100,xo(2,),,Tc)
gtgt Tc 310 gtgt t,x ode45(df,0
10,xo(3,),,Tc)
gtgt axplotyy(t,x(,1),t,x(,2)) gtgt xlabel('Time
h') gtgt ylabel(ax(2),'Temperature K') gtgt
ylabel(ax(1),'Concentration kmol/m3')
8
Matlab Exercise cstr.m
function f cstr(x,Tc) ko 3.493e7 E
11843 H -5960 rhoCp 500 UA 150 R
1.987 V 1 q 1 Caf 10 Ti 298 Ca
x(1) T x(2) f(1) q/V(Caf - Ca) -
koexp(-E/R/T)Ca f(2) q/V(Ti - T)
-H/rhoCpkoexp(-E/R/T)Ca UA/rhoCp/V(Tc-T) f
f'