Using Matlab Solving Ordinary Differential Equations Analytically (by T. Gu)

Matlab can solve linear ODE's or ODE systems with or without initial/boundary conditions. Do not expect Matlab to solve nonlinear ODE's which typically have no analytical solutions. Matlab version R12 and later treats y as default function symbol and t as default variable symbol. Earlier version treats x as default symbol.

Example 1.  d2y/dx2 -2dy/dx -3y=x2
>> dsolve('D2y - 2*Dy - 3*y=x^2', 'x')
ans =
-14/27+4/9*x-1/3*x^2+C1*exp(3*x)+C2*exp(-x)
Example 2. d2y/dx2 -2dy/dx -3y=x2, with y(0)=0, and y(1)=1. 
>> dsolve('D2y - 2*Dy - 3*y=x^2','y(0)=0, y(1)=1','x')
ans =
-14/27+4/9*x-1/3*x^2+2/27*(7*exp(-1)-19)/(exp(-1)-exp(3))*exp(3*x)-2/27*(-19+7*exp(3))/(exp(-1)-exp(3))*exp(-x)

Example 3. d2y/dx2 -2dy/dx -3y=x2, with y(0)=0, and dy/dx =1 at x=1
>> dsolve('D2y - 2*Dy - 3*y=x^2','y(0)=0, Dy(1)=1','x')
ans =
-1/3*x^2+4/9*x-14/27+1/9*(-11+14*exp(3))/(3*exp(3)+exp(-1))*exp(-x)
+1/27*(33+1 4*exp(-1))/(3*exp(3)+exp(-1))*exp(3*x)

Example 4.  d2y/dx2 -2dy/dx -3y=0, with y(0)=a, and y(1)=b.
>> dsolve('D2y-2*Dy-3*y=0','y(0)=a, y(1)=b')
ans =
(exp(3)*a-b)/(-exp(-1)+exp(3))*exp(-t)-(-b+a*exp(-1))/(-exp(-1)+exp(3))*exp(3*t)
>> pretty(ans)
              (exp(3) a - b) exp(-t)   (-b + a exp(-1)) exp(3 t)
              ---------------------- - -------------------------
                -exp(-1) + exp(3)          -exp(-1) + exp(3)

Example 5.  d2y/dt2+y=4cos(t), y(pi/2)=2pi, dy/dt=-3 at t=pi/2
>> y=dsolve('D2y+y=4*cos(t)' , 'y(pi/2)=2*pi, Dy(pi/2)=-3')
y =
-2*sin(t)^2*cos(t)+(2*sin(t)*cos(t)+2*t)*sin(t)+5*cos(t)+pi*sin(t)
>> simplify(y)
ans =
2*sin(t)*t+5*cos(t)+pi*sin(t)

Example 6.  Solving two ODE's with two initial conditions:
df/dx =3f+4g, dg/dx =-4f+3g, with f(0)=1, g(0)=1
>> [f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','f(0)=0,g(0)=1','x')
f =
exp(3*x)*sin(4*x)
g =
exp(3*x)*cos(4*x)

Example 7.  Bessel Equation: x2(d2y/dx2) +x(dy/dx) + (x2-v2)y=0
                 Solution: y=c1Jv(x) + c2Yv(x)
>> y=dsolve('x^2*D2y+x*Dy+(x^2 - v^2)*y =0','x')
y =
C1*besselj(v,x)+C2*bessely(v,x)


Using Matlab Solving Algebraic Equations Analytically
Matlab can be used to solve nonlinear algebraic equations or equation systems. Sometimes only numeric answers are returned. Sometimes, Matlab fails.
Example 1:
>> solve('-x^3 +3*x =0')
ans =
[       0]
[-3^(1/2)]
[ 3^(1/2)]

>> numeric(ans)
ans =
         0
   -1.7321
    1.7321

Or, use the command "roots" for polynomials
>> roots([-1,0,3,0])
ans =
         0
    1.7321
   -1.7321


Example 2:
>> solve('ln(x)+a^2=b')
ans =
exp(-a^2+b)

Example 3:
>> T=solve('ln(T)+a^2=b','T')
T =
exp(-a^2+b)

Example 4: Solving 2 equations simultaneously (Output from version 5. Version R12 fails to solve it!)
>> e1='x1^2 + x2 =10';
>> e2='x1-x2=2';
>> solve(e1,e2,'x1,x2')
ans =
x1 = -4, x2 = -6
x2 = 1, x1 = 3  









Using MATLAB to Solve for Friction Factor Using Colebrook Formula function y=cole(f) eoverD=0.000375; Re=1.37e+4; y=1/sqrt(f)+2*log10(eoverD/3.7 + 2.51/Re/sqrt(f)); Save the text above as cole.m in c:\temp In Matlab, type cd c:\temp f=fsolve('cole', 0.02) The following will show up on screen » cd c:\temp » f=fsolve('cole', 0.02) f =     0.0291 Comment:  Using f=solve('...') does not work. Need to use fsolve with an initial guess. Using fsolve for a nonlinear equation system The problem below is from a multiple pipe system in ChE345 (Fluid Mechanics.) function z=mflow(v) z(1)=v(1)-v(2)-v(3); z(2)=322 -v(1)^2-0.4*v(3)^2; z(3)=258 -v(1)^2-0.5*v(2)^2; Save the text above as a file called mflow.m in c:\temp In Matlab, type the following three command lines: cd c:\temp x0=[1,1,1] fsolve('mflow',x0) The screen shows >> cd c:\temp >> x0=[1,1,1] x0 =      1     1     1 >> fsolve('mflow',x0) ans =    15.9327    2.8802   13.0526 The same problem can be solved using the solve command: » e1='x1 = x2 + x3'; » e2='322 = x1^2+0.4*x3^2'; » e3='258 = x1^2+0.5*x2^2'; » solve(e1,e2,e3,'x1,x2,x3') ans = x1 = -5.665593102967009, x3 = -26.92123022781521, x2 = 21.25563712484820  x1 = 5.665593102967009, x3 = 26.92123022781521, x2 = -21.25563712484820   x1 = 15.93274220916813, x3 = 13.05255968155214, x2 = 2.880182527615988    x1 = -15.93274220916813, x3 = -13.05255968155214, x2 = -2.880182527615988 "solve" is not as reliable as "fsolve" in MATLAB. Another example using fsolve to solve 6 nonlinear equations This is for Problem 8.64 on p. 374 of ChE345 (Fluid Mechanics) textbook. How to solve Homework Problem 8.64* using Matlab's fsolve? 1. Create a function m-file named hw864.m with the follow lines in it function z=fcn(x) z(1)=x(1)-0.64*x(2)-0.64*x(3);    %x(1)=v1, x(2)=v2, x(3)=v3 z(2)=40-102*x(4)*x(1)^2 - 127*x(5)*x(2)^2; %x(4)=f1, x(5)=f2, x(6)=f3 z(3)=60-102*x(4)*x(1)^2 - 255*x(6)*x(3)^2; z(4)=1/sqrt(x(4))+2.0*log10(1.22e-4 +2.81e-5/x(1)/sqrt(x(4)));   %Colebrook for f1 z(5)=1/sqrt(x(5))+2.0*log10(1.52e-4 +3.52e-5/x(2)/sqrt(x(5)));   %Colebrook for f2 z(6)=1/sqrt(x(6))+2.0*log10(1.52e-4 +3.52e-5/x(3)/sqrt(x(6)));   %Colebrook for f3   2. Save hw864.m in c:\temp 3. Do the following in Matlab: >> cd c:\temp >> x0=[1,1,1,0.02,0.02,0.02]; >> fsolve('hw864',x0) ans =     3.5001    2.6915    2.7774    0.0179    0.0192    0.0191 The answers are: v1=3.500 m/s, v2=2.692 m/s, v3=2.777 m/s f1=0.0179, f2=0.0192, f3=0.0191
Using Matlab for Laplace Transform
(The RCENT version of MATLAB V6.0.0.88/Release 12 forgot to place a copy of laplace.m in its toolbox\matlab\specfun directory. Copy it from its toolbox\symbolic\@sym directory. Otherwise, laplace command is invalid. Same for the ilaplace command. Version 5 uses the name invlaplace.)
Example 1.
>> laplace('exp(-a*t)*t')
ans =
1/(s+a)^2
>> ilaplace('1/(s+a)^2')            %verify the previous result
ans =
exp(-a*t)*t

Example 2.
>> laplace('exp(-a*t)*cos(b*t)')
ans =
(s+a)/((s+a)^2+b^2)
EDU» pretty             
 
                                     s + a
                                 -------------
                                        2    2
                                 (s + a)  + b
 
>> ilaplace('(s+a)/((s+a)^2 + b^2)')
ans =
exp(-a*t)*cos(b*t)


Differentiation and Integration Using Matlab:
Note: Integration (without limits) results omit the integration constant in the answer.
Example 1.
>> diff('x^3 + ln(x)*sin(x) + b/x')
ans =
3*x^2+1/x*sin(x)+log(x)*cos(x)-b/x^2
>> int('3*x^2+1/x*sin(x)+log(x)*cos(x)-b/x^2')         %verify
ans =
log(x)*sin(x)+(x^4+b)/x

Example 2. Partial Derivative: Take partial derivative with respect to t.
>> diff('x*t^3+ln(t)','t')   %specify t as the variable. x will be treated as a const.
ans =
3*x*t^2+1/t
>> int('3*x*t^2+1/t','t')         %verify
ans =
x*t^3+log(t)

Example 3. Integrate ln(x)+1/(x+1) from x=1 to x=2. 
>> int('ln(x) + 1/(x+1)', 1, 2)
ans =
-1+log(2)+log(3)
>> eval(ans)
ans =
    0.7918