gu_examp

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

Click here to see numerical tools with graphical user interface developed 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.  d²y/dx² -2dy/dx -3y=x²>> 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. d²y/dx² -2dy/dx -3y=x², 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. d²y/dx² -2dy/dx -3y=x², 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.  d²y/dx² -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: x²(d²y/dx²) +x(dy/dx) + (x²-v²)y=0

                 Solution: y=c₁J_v(x) + c₂Y_v(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
>> e1='x1^2 + x2 =10';
>> e2='x1-x2=2';
>> [x1,x2]=solve(e1,e2)
x1 =
[ -4]
[  3]
x2 =
[ -6]
[  1]

You may also use a single command line shown below to get the same answer.
>> [x1,x2]=solve('x1^2+x2=10','x1-x2=2')

Using MATLAB to Solve for Friction Factor Using Colebrook Formula
Save the following 6 lines into friction.m m-file.

function dummy %this allows you to pack the 5 lines below into a single m-file.
fsolve(@cole, 0.02) %fsolve looks for cole function below. If missing, it looks in current directory.
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));

Run friction.m in Matlab
>>friction
f =
    0.0291


Using fsolve for a nonlinear equation system

The problem below is from a multiple pipe system in ChE345 (Fluid Mechanics.)
Save the 7 lines below as a file named fluids.m to be run in Matlab.

function dummy  %This allows you to pack the 6 lines below into a single m-file
x0=[1,1,1];
fsolve(@mflow,x0)  %fsolve looks for mflow function below. If missing, it looks in current directory.
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;

Run the m-file to get three velocities. 
>> fluids
ans =
   15.9327    2.8802   13.0526

The same problem can be solved using the solve command.
Save the 4 lines in equations.m m-file.
e1='x1 = x2 + x3';
e2='322 = x1^2+0.4*x3^2';
e3='258 = x1^2+0.5*x2^2';
[x1,x2,x3]=solve(e1,e2,e3)

>>equations
x1 =
 -15.932742209168131268551587928289
 -5.6655931029670086917448152900978
  5.6655931029670086917448152900978
  15.932742209168131268551587928289
 
x2 =
 -2.8801825276159898268394750195246
  21.255637124848206348954897950637
 -21.255637124848206348954897950637
  2.8801825276159898268394750195246
  
x3 =
 -13.052559681552141441712112908765
 -26.921230227815215040699713240735
  26.921230227815215040699713240735
  13.052559681552141441712112908765

"solve" is not as reliable as "fsolve" in MATLAB.
When you copy text from WORD or html to Matlab, be careful! You may carry over some hidden symbols. Use a notepad file to drop them. 


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?

Save the 10 text lines below as solution.m file.

function dummy
x0=[1,1,1,0.02,0.02,0.02];
fsolve(@fcn,x0)
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

Run solution.m in Mablat

>>solution
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

Factoring a function or number
>> factor(12)
ans =
2 2 3
>> syms x y % declare that x and y are symbols
>> factor(x^2-y^2)
ans =
(x-y)*(x+y)
>> factor(x^2-1)
ans =
(x-1)*(x+1)

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)
>> pretty             

                                     s + a

                                 -------------

                                        2    2

                                 (s + a)  + b

 >> ilaplace('(s+a)/((s+a)^2 + b^2)')
ans =
exp(-a*t)*cos(b*t)

Example 3.
>> sym('n','positive')            %This line is needed!!!
ans =
n
>> laplace('t^n')
ans =
s^(-n-1)*gamma(n+1)
>> ilaplace('s^(-n-1)*gamma(n+1)')
ans =
t^n

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