Matlab符号计算与⽅程组求解
⼀、符号计算
1、符号计算特点
1、计算精确:符号计算基于数学公式、定理并通过⼀系列推理、演绎得到⽅程的解或者数学表达式的值。对操作对象不进⾏离散化和近似化处理。
2、可应⽤范围有限:实际科研和⽣产中遇到的问题绝⼤多数都⽆法获得精确的符号解,这时我们不得不求助数值计算。
3、对待符号计算态度:⽤其来完成公式推导和解决简单的对计算时效性要求不⾼的问题,综合符号计算和数值计算各⾃的优点,视问题特点混合使⽤符号计算和数值计算。
2、符号对象和符号表达式
1、符号对象的创建:要⽣成⼀个符号对象,可以利⽤sym以及syms函数,sym可以⽣成单个符号对象,⽽syms可以⽣成多个符号对象,符号对象的运算是完全精确的,没有舍⼊误差。例如:a = sym(‘5’);  syms b c d;。
2、符号表达式:创建了符号对象,我们就可以创建各种各样的符号表达式。譬如,创建符号变量a, b, c后如下都是符号表达式:
z1 = a+b+c;
z2 = sin(a+b+c);
z3= a^b*gamma(c);
确定⼀个符号表达式中的符号变量可以⽤findsym函数:
findsym(expr);
findsym(expr, n);
第⼀种⽤法是确认表达式expr中所有⾃由符号变量,第⼆种⽤法是从表达式expr中确认出距离x最近的n个符号变量。这个最近距离指的是变量第⼀个字符和x的ASCII码值之差的绝对值,差绝对值相同时,ASCII码值⼤的字符优先。findsym是⽼版本中的⽅法,新版本将使⽤symvar函数代替:symvar(expr);。
3、运算符:MATLAB采⽤了重载(Overload)技术,使得⽤来构成符号表达式的运算符,⽆论在拼写还是在使⽤⽅法上,都与数值计算中的算符完全相同。譬如“+”,“-”,“*”,“\”,“/”,“^”等。
符号对象的⽐较中,没有“⼤于”、“⼤于等于”,“⼩于”,“⼩于等于”的概念,⽽只有是否“等于”的概念,
即“==”与“~=”。如果要判断两个符号数值的⼤⼩⼀般来说有两种办法,⼀种是利⽤double将其转化成数值型的,另⼀种是利⽤
sort+“==”或“~=”。譬如下⾯代码:
1. >> a = sym('2');
2. >> b= sym('3');
3. >>double(a)<double(b)
4. ans =
5. 1
6. >> sa = sort([b,a])
7. sa =
8. [ 2, 3]
9. >>a==sa(1)
0. ans =
1. 1
从上述代码可以看出,上述两种⽅法都间接实现了判断⼤⼩。
4、符号计算与数值计算结合:利⽤符号计算得到结果时,有时需要将其转化成数值型的以便后续数值计算利⽤。通过符号计算得到⼀个表达式时,想把它转化成关于其中某个变量的数值函数。
很多时候我们需要求符号表达式在不同的参数值下的具体值,说通俗点就是如何把具体的参数代⼊符号表达式。这时候可以利⽤eval 和subs函数或者转化成匿名函数。
程序实例:
1. clc;
2. clear;
3.
4. % 将符号变量转化成数值
5. format long;
6. % 可供使⽤的⼀些⽅法
7. % 这⾥a是符号变量,a1--a4是数值变量
8. a = vpa(pi,30); %3.14159265358979323846264338328
9. a1 = double(a); %3.141592653589793
0. a2 = eval(a); %3.141592653589793
1. a3 = single(a); %3.1415927
2. a4 = int8(a); % 3
13.
14.
5. % 很多时候我们需要求符号表达式在不同的参数值下的具体值,
6. % 说通俗点就是如何把具体的参数代⼊符号表达式。这时候可以利⽤eval和subs函数或者转化成匿名函数
17.
8. % 如下的例⼦:求函数:f =sin( x^x / (x^2/exp(x)) ); 其⼆阶导数在x = 1处的值。
9. syms x;
0. f = sin(x^x / x^2/exp(x));
1. % 利⽤符号计算求f(x)的⼆阶导数
2. % diff函数⽤于求导数或者向量和矩阵的⽐较。
3. % 如果输⼊⼀个长度为n的⼀维向量,则该函数将会返回长度为n-1的向量,向量的值是原向量相邻元素的差
4. d2f = diff(f, x, 2);
25.
6. % 第⼀种⽅法:利⽤subs函数求d2f在x=1时的值。
7. d2fx1 = subs(d2f, x, 1); % d2fx1 = 2.2082
28.
9. % 第⼆种⽅法:x赋值1后,利⽤eval函数求d2f在x = 1时的值
0. x = 1;
1. d2_fx1 = eval(d2f); %d2_fx1 =
2.2082
32.
3. % 第三种⽅法:将d2f转化成匿名函数,求其在x = 1时的值
4. % vectorize的含义就是将乘转成点乘等。 '*' -> '.*'; '/' -> './'; '^' -> '.^';最后再将替换结果中的“..”删除⼀个"."。
5. F = eval(['@(x)',vectorize(char(d2f))]);
6. F(1) % ans= 2.2082
⼆、极限、导数和级数的符号计算
1、求的极限。
1. clc;
2. clear;
3.
4. % 求极限
5. syms n;
6. % limit函数⽤于求极限运算。假如y=f(x),limit(y)表⽰x→0时的极限,limit(y,x,a)表⽰x→a时的极限
7. % MATLAB中gamma函数即数学上的gamma函数,当n为正整数时,有如下性质:gamma(n+1) = n!
8. limit( n^(n+1/2) /( exp(n)*gamma(n+1)), n,inf) % ans= 1/(2*pi)^(1/2)
2、求的导数,。
1. % 求导数
2. syms a t x;
3. f = [a, t*log(x); sqrt(t), x^2+3*x];
4. % 矩阵f对t的⼀阶导数
5. dfdt = diff(f,t); % dfdt = [ 0, log(x); 1/(2*t^(1/2)), 0 ];
6.
7. % 矩阵f对x的⼆阶导数,由于是x,⽽f中含有x变量,故x可以省略
8. dfdx2 = diff(f,2); % dfdx2 =[ 0, -t/x^2; 0, 2 ];
9.
0. %求⼆阶混合导数
1. dfdtdx = diff(diff(f,t),x); %dfdtdx =[ 0, 1/x; 0, 0 ];
3、求⽆穷级数:
% 求级数
syms k;
f1 = symsum( (k-2)/2^k, k, 3, inf);  % f1 = 1/2
A = [1/(2*k+1)^2, (-1)^k/3^k];
f2 = symsum(A, k, 1, inf); % f2= [ pi^2/8 - 1, -1/4]
三、求解⽅程
求解⽅程通常有两种⽅法,符号求解和数值求解。
1、solve
通常在不确定⽅程是否有符号解的时候,推荐先使⽤solve进⾏尝试,因为solve相⽐于数值求解来说,它不需要提供初值,并且⼀般情况下能够得到⽅程的所有解。对于⼀些简单的超越⽅程,solve还能够⾃动调⽤数值计算系统给出⼀个数值解。
solve的调⽤形式:
sol = solve(eq)
sol = solve(eq, var)
sol = solve(eq1, eq2, …, eqn)
sol = solve(eq1, eq2, …, eqn, var1, var2, …, varn)
eq为符号表达式,var为指定的要求解的变量。如果不声明要求解的变量(第⼀和第三种形式),则matlab⾃动按默认变量进⾏求解,默认变量可以由symvar(eq)确定。
例:求解⽅程组: x+y=1, x-11y=5
1. % 使⽤solve函数求解⽅程
2. % 例:求解⽅程组: x+y=1, x-11y=5
3. clc;
4. clear;
5.
6. % 声明符号变量
diff函数7. syms x y;
8. % 定义⽅程
9. eq1 = x + y -1;
0. eq2 = x - 11*y - 5;
1. % 调⽤solve求解⽅程组
2. % (solve函数的参数包括⽅程表达式,以及要求解的变量,这⾥变量是可选参数,不指定时matlab⾃动按默认变量进⾏求解)
3. % 这时候solve求得的解通过结构体的形式赋值给sol,然后再通过x=sol.x和y=sol.y分别赋值给x和y。
4. sol = solve(eq1, eq2, x, y);
5. x = sol.x;
6. y = sol.y;
7. % 也可以直接使⽤:
8. % [x,y] = solve(eq1, eq2, x, y);
19.
0. % 使⽤double将符号解转换为数值解
1. value_x = double(x);
2. value_y = double(y);