Wednesday, 28 November 2007

Symbolic Operations in Matlab

Background Note:

Up to this point we have done purely numerical operations in Matlab. This is the traditional approach for solving problems on the computer, and much of the current computational work done in industry follows this approach. However, a relatively new capability (last 5-10 years or so) that is gaining popularity involves the use of symbolic operations directly on the computer. Computer programs that use these symbolic techniques are often referred to as computer algebra systems, and they are becoming quite powerful. In particular, Matlab’s Symbolic Math Toolbox offers this general capability within the Matlab environment by providing an interface to the Maple code (Maple is a commercial software package that emphasizes symbolic manipulations or computer algebra techniques). This set of labs will introduce some of the basic symbolic capability in Matlab, eventually leading up to the analytical solution of Ordinary Differential Equations directly on the computer. Many of the exercises given here are derived from the Matlab’s User’s Guide and you are encouraged to refer to that manual for further examples and guidance.

Define Some Symbolic Variables:

First define a bunch of symbolic variables: syms a1 b1 c1 a2 b2 c2 d2 A B C1 C2 x y t

Now let’s form some symbolic functions:

Mathematical Function
Matlab Syntax for Function

f1 = a1 + b1*x + c1*x^2

f2 = a2 + b2*x + c2*x^2 + d2*x^3

g = exp(A*t)*(C1*cos(B*t)+C2*sin(B*t))

u = 2*x*y^2 + sin(x+y)

Differentiation of Symbolic Functions:

Differentiate each of the above functions using Matlab’s diff command. For example, try diff(f1). What happened? Is the result correct? Do the same thing for the other functions. What happens when you type diff(u)? How do you get Matlab to take the partial derivative with respect to the variable y? Type help sym/diff to see the various options associated with taking symbolic derivatives.

We can also take higher order derivatives. For example diff(f1,2) takes the second derivative of f1(x). Try this for all the functions!

Can you show that the mixed partial derivatives of u(x,y) are equal? Take du/dx and du/dy as above and save the results as string variables. For example, try typing

dudx = diff(u,x); dudy = diff(u,y);.

Now take the first derivative of these symbolic variables, as follows:

uyx = diff(dudx,y); uxy = diff(dudy,x);

Are these latter two variables the same? They should be. Pretty neat stuff, don’t you think!!!

Integration of Symbolic Functions:

Using Matlab’s int function, perform the following integrals using f1(x):

Mathematical Operation
Matlab Syntax for Operation

I1 = int(f1)

I2 = int(f1, 0, 5)

Do these make sense? Check out the subs command. Now let a1 = 1, b1 = -2, and c1 = 1, and substitute these numerical constants into the symbolic expressions for I1 and I2, as follows:.

I1 = subs(I1,{a1 b1 c1},{1 -2 1}), I2 = subs(I2,{a1 b1 c1},{1 -2 1})

What are the values for the integrals I1 and I2? Why is I1 a function of x and I2 is simply a scalar number? Note also that the curly brackets in Matlab refer to cell arrays. The sequence given by the set, {a1 b1 c1},{1 -2 1}, replaces the numerical variables into the symbolic variables in sequence. The cell array simply let’s us make all three substitutions in one call to the subs command.

Try integrating and substituting values for the constants in the other functions. Note that with u(x,y), you can only integrate with respect to one variable at a time. For example,

Mathematical Operation
Matlab Syntax for Operation

Iu = int(int(u,x,0,pi),y,0,2*pi)

Pretty impressive!!! Could you do this integral by hand this fast? I certainly can’t…

Plotting Symbolic Functions:

Certainly we can also plot symbolic relationships. Since we have already defined the constants for use with f1(x), let’s use this function to do some simple plots. In particular, let’s plot f1(x), df1/dx, and over the range . We can do this as follows:

F1 = subs(f1,{a1 b1 c1},{1 -2 1}); D1 = diff(F1); I1 = int(F1);
Nx = 51; xp = linspace(0,5,Nx);
F1p = double(subs(F1,x,xp));
D1p = double(subs(D1,x,xp));
I1p = double(subs(I1,x,xp));
title(‘Evaluating and Plotting Symbolic Functions’)
xlabel(‘X Values’),ylabel(‘Various Function Values’)
legend(‘Function’,’First Derivative’,’Integral’)

Solving Algebraic Equations:

We can also solve algebraic equations with Matlab. For example, what if we wanted to know the value of x where f1(x) = 3? With the above constants this problem becomes: find x such that

In Matlab, we can use the solve command: solve(x^2-2*x-2) or solve(F1 - 3).

We can also solve systems of algebraic equations. For example, the analytical solution to the following 2x2 system:

is given by:

S = solve(3*x + 2*y - 3,-2*x + y + 7); % this solves a simple set of 2x2 equations

S.x, S.y % this prints x and y to the screen

Solving Ordinary Differential Equations:

The real goal of the above discussion and examples with symbolic variables is to give enough background so that we can solve ODEs analytically on the computer. This will be pretty powerful capability if we can actually do this with a set of relatively simple commands. As a test, let’s give Matlab’s dsolve command a workout for a few different types of ODEs that we have treated thus far in the semester (be sure to type help dsolve to get a good idea of the various forms that can be used).

a. Solve the first order IVP:

dsolve(‘Dy + y/x = (2/x)*exp(2*x)’,’x’) % this gives the general solution

dsolve(‘Dy + y/x = (2/x)*exp(2*x)’,’y(1) = 0’,’x’) % this gives the unique solution

b. Solve the second order IVP:

Sh = dsolve('D2y - 4*y = 0','x'), pretty(simple(Sh)) % homogeneous soln
Sg = dsolve('D2y - 4*y = 4*x^2','x'), pretty(simple(Sg)) % general soln
Su = dsolve('D2y - 4*y = 4*x^2','y(0) = -1/2','Dy(0) = 4','x') % unique soln

c. Solve the second order IVP:

Q1 = dsolve(‘D2y + 6*Dy +13*y = 10*sin(5*t)’, ’y(0) = 0’,’Dy(0) = 0’,’t’)

d. Re-solve the problem in Part c as a system of two 1st order ODEs, where z1 = y and z2 = dy/dt:

Q2 = dsolve(‘Dz1 =z2’,’Dz2 = -13*z1 -6*z2 + 10*sin(5*t)’, ’z1(0) = 0’,’z2(0) = 0’,’t’)
pretty(simple(Q2.z1)), pretty(simple(Q2.z2))

Final Note:

There is a lot of good stuff here. In particular, there are several examples that should make your life much easier in this course and in several of your other technical classes. You should do your best to understand the basic capability that is illustrated here -- it represents some pretty powerful mathematical analysis capability. You should also note that some of the examples from earlier in the semester (the Two Salty Tanks problem for example) give other illustrations of the use of computer algebra to solve a coupled set of first order differential equations for a real problem of interest. Also, a final Matlab demo will be given at the end of the semester illustrating how to take Laplace transforms and inverse Laplace transforms analytically using Matlab’s symbolic processing capability.


No comments:

My photo
London, United Kingdom

Facebook & Twitter