一.符号函数
暂时只介绍代码部分原理暂不解释。
clear;clc
%% 符号方程的创建
syms a x
y = a*x+x^2
%% 化简
syms a
y=(cot(a/2)-tan(a/2))*(1+tan(a)*tan(a/2))
simplify(y) % 化简
%% 因式分解
factor(12) % 对常数进行因式分解
% 2 2 3
syms m n x
y = -24*m^2*x-16*n^2*x
factor(y) % 因式分解
% [ -8, x, 3*m^2 + 2*n^2]
y1=m^3-n^3
factor(y1)
% [ m - n, m^2 + m*n + n^2]
%% 多项式展开
syms a x
y = a*(x^2-a)^2+(x-2)
expand(y) % 多项式展开
% a^3 - 2*a^2*x^2 + a*x^4 + x - 2
%% 合并
syms x y
z = (x+y)^2*y+5*y*x-2*x^3
% expand(z) - 2*x^3 + x^2*y + 2*x*y^2 + 5*x*y + y^3
collect(z,x)
% y*x^2 - 2*x^3 + (2*y^2 + 5*y)*x + y^3
collect(z,y)
% y^3 + 2*x*y^2 + (x^2 + 5*x)*y - 2*x^3
%% 计算分子与分母
% [z1,z2] = numden(2.5) % 会报错,因为numden的输入变量不能是数值,只能是符号变量
% ans = sym(2.5); % sym函数可以将数值2.5转换为符号
[z1,z2] = numden(sym(2.5)) % 对常数计算分子与分母
syms x y
z = 1/x*y+x/(x^2-2*y)
[z1,z2] = numden(z) %z1分子,z2分母
%% 让结果显示的更加自然
syms x y
M = (1/x*y+x/(x^2-2*y)-x^2/(3+y)^2)^2;
expand(M)
2.符号函数求导(diff函数用法)
% 一元函数的导数
syms x
y = x^4-5*x^2+6
diff(y) %求一阶导数
diff(y,2) %求二阶导数
% 多元函数的导数
syms x1 x2 x3
y1 = x1^5*x2+x2*x3-x1^2*x3
py1 = diff(y1,x1,1) % 对x1求一阶偏导
% 5*x2*x1^4 - 2*x3*x1
py2 = diff(y1,x1,2) % 对x1求二阶偏导
% 20*x2*x1^3 - 2*x3
py3 = diff(y1,x1,x2) % 先对x1求偏导,再对x2求偏导
% 5*x1^4
py4 = diff(y1,x2,x1) % 先对x2求偏导,再对x1求偏导
% 5*x1^4
%% 注意,如果diff函数作用的对象不是符号函数,而是矩阵,那么对应的功能是求差分。
A=[4 5 6 3 2 1];
diff(A) % 求向量A的一阶差分 1 1 -3 -1 -1
diff(A,2) % 在一阶差分的基础上再差分一次 0 -4 2 0
A=[4 5 6;
7 4 2;
5 6 2]
A1=diff(A) % 下一行减去上一行求一阶差分
% 3 -1 -4
% -2 2 0
A2=diff(A,2) % 下一行减去上一行求二阶差分(在一阶差分的基础上再差分一次)
% -5 3 4
A3=diff(A,2,1) % 最后面的1表示在行上进行差分(在列的方向上进行差分)
% -5 3 4
A4=diff(A,1,2) % 后一列减去前一列求一阶差分, 最后面的2表示在列上进行差分(在行的方向上进行差分)
% 1 1
% -3 -2
% 1 -4
A4=diff(A,2,2) % 后一列减去前一列求二阶差分
% 0
% 1
% -5
二.不定积分
%% 不定积分
clear;clc
syms x
y = x^2
int(y,x)
% x^3/3 注意,Matlab计算时不会给我们加上常数C
syms x
y = 1/x
int(y,x)
% log(x) 注意,Matlab计算1/x形式的不定积分时不会给我们加上绝对值~
syms x
y = x^2 / (1+x^2)
int(y,x)
% x - atan(x)
syms x
y = 1/(exp(x)+1)
int(y,x)
% x - log(exp(x) + 1)
syms x a
y = 1/sqrt(x^2-a^2)
int(y,x)
% log(x + (x^2 - a^2)^(1/2))
%% 定积分
syms x
y = sin(x)
int(y,x,0,pi/2)
% 1
syms x a b
y = exp(x)
int(y,x,a,b)
% exp(b) - exp(a)
syms x
y = (sin(x))^2 / x^2
b=int(y,x,0,+inf)
% pi/2
% 注意,不是所有的函数都可以利用int函数计算出最后的结果,例如:
syms x
y = 1 / exp(x) * log(x+2*x^2+sin(x))
int(y,x,0,4)
% int(exp(-x)*log(x + sin(x) + 2*x^2), x, 0, 4)
% 我们可以计算数值积分:数值积分可用于求定积分的近似值。在数值分析中,数值积分是计算定积分数值的方法和理论。
% 在数学分析中,给定函数的定积分的计算不总是可行的,许多定积分不能用已知的积分公式得到精确值。
y = @(x) 1 ./ exp(x) .* log(x+2.*x.^2+sin(x)) % 注意,写成函数句柄时,要用点乘或者点除
integral(y,0,4)
xx = 0:0.1:4;
yy = 1 ./ exp(xx) .* log(xx+2*xx.^2+sin(xx));
plot(xx,yy,'-')
三.方程解和方程组解
%% solve函数
%% 例题1: 求解单变量方程
clear;clc
syms x
answ = solve(sin(x) == 1, x) % 注意:这里的等号一定要有两个,一个等号表示赋值,两个等号才表示左右两边相等
answ = solve(sin(x) == 1) % 只有一个符号变量x,所以可以不指定未知数
% 也可以这样写
syms x
eqn = (sin(x) == 1); % eqn = sin(x) == 1;
answ = solve(eqn, x)
% 因为三角函数是周期函数,如果要得到所有的解,则需要加上条件
[answ, params, condions] = solve(eqn, x, 'ReturnConditions', true)
%% 例题2: 多变量方程求解
syms a b c x
eqn = (a*x^2 + b*x + c == 0);
answ1 = solve(eqn, x) % 将x视为未知数求解
answ2 = solve(eqn, a) % 将a视为未知数求解
%% 例题3:方程组求解
syms u v a
eqn = [2*u + v == a, u - v == 1];
answ = solve(eqn, [u, v])
answ.u
answ.v
[answ_u, answ_v] = solve(eqn, [u, v])
%% solve 可能会警告
syms x
eqn = (sin(x) == x^2 - 1);
solve(eqn, x) % 警告: Cannot solve symbolically. Returning a numeric approximation instead.
% 画图看看
fplot(sin(x), [-2 2]) % fplot函数可绘制表达式的图形
hold on
fplot(x^2 - 1, [-2 2])
%% vpasolve函数求解
% 用vpasolve函数指定求[0 2]上的解
syms x
eqn = sin(x) == x^2 - 1;
vpasolve(eqn, x, [0 2]) %只会返回一个解
vpasolve(eqn, x, 'random', true) %会返回多个解
vpasolve(eqn, x, -5) % 给定搜索的起始点
%% 来看一个更复杂的例子
syms x y
eqn = [x^2 - 2*x - 3*x*y == 10, y^4 == exp(-2*x/3*y)]
[answ_x, answ_y] = vpasolve(eqn, [x, y], 'random', true)
%绘制隐函数图像
fimplicit(x^2 - 2*x - 3*x*y == 10, [-10 10],'r') % R2016b版本之后才有
hold on
fimplicit(y^4 == exp(-2*x/3*y), [-10 10],'b') % R2016b版本之后才有
[answ_x, answ_y] = vpasolve(eqn, [x, y],[-4 -1;1 5]) % 指定搜索的范围:x位于[-4 -1], y位于[1 5]
hold on
% plot(double(answ_x), double(answ_y),'ko', 'MarkerSize',10) % double可以将我们的符号变量转换为数值变量
%% fsolve函数(求解功能最为强大哦)
% fsolve是Matlab优化工具箱中的一个函数,可专门用来求解特别复杂的方程和方程组
x0 = [0,0]; % 初始值
result_x = fsolve(@my_fun,x0)
% 当然你也可以用vpasolve函数试试
clear; clc
syms x1 x2
eqn = [exp(-exp(-(x1+x2))) - x2*(1+x1^2) == 0, x1*cos(x2) + x2*sin(x1) - 0.5 == 0]
[answ_x1, answ_x2] = vpasolve(eqn, [x1, x2], [0 0])
function F = my_fun(x)
F(1) = exp(-exp(-(x(1)+x(2)))) - x(2)*(1+x(1)^2);
F(2) = x(1)*cos(x(2)) + x(2)*sin(x(1)) - 0.5;
end
四.微分方程解
1.解析解
clc;clear
%%
syms y(x) a
eqn = (y - diff(y,x) == a*x);
dsolve(eqn)
syms y(x)
eqn = (y - diff(y,x) == 2*x);
cond = (y(0) == 3);
dsolve(eqn,cond)
syms y(x)
eqn = (diff(y,x,2) + 4 *diff(y,x) + 29*y == 0);
Dy = diff(y,x); % 定义变量Dy为y的一阶导数
cond = [(y(0) == 0) ,(Dy(0) ==15)] ; % 有两个条件,可以写到一个向量中保存
dsolve(eqn,cond)
%% 微分方程组
syms x(t) y(t) z(t)
eqn1 = (diff(x,t) == 2*x-3*y+3*z+t);
eqn2 = (diff(y,t) == 4*x-5*y+3*z+t);
eqn3 = (diff(z,t) == 4*x-4*y+2*z+t);
eqns = [eqn1 eqn2 eqn3];
[x,y,z] = dsolve(eqns)
2.数值解
(1)一阶微分方程
%% 调用ode45函数求解微分方程df1,自变量为x,范围为[0,2], 初始值y(0)=3 ; 因变量为y
clear;clc
[x,y] = ode45('df1',[0,2],3); % [x,y] = ode45(@df1,[0,2],3);
%% 设定相对误差和绝对误差,这样可以提高微分方程数值解的精度
options = odeset('reltol',1e-4,'abstol',1e-8);
[x,y] = ode45('df1',[0,2],3,options);
%% 如果觉得x的间隔不够小,我们可以指定要求解的位置
[x,y] = ode45('df1',[0:0.001:2],3,options);
function dy = df1(x,y)
% 微分方程:y-y'=2x(函数名称可以任意取)
dy = y - 2*x; % 写成标准形式 y' = y - 2x
% 注意函数的返回值一定是因变量y的一阶导数
% 函数的输入有两个,分别是自变量x和因变量y
end
(2)一阶微分方程组
%% 在真正比赛中,我们往往会倾向于先去看有无解析解,如果没有解析解再考虑数值解
% [y1 y2 y3] = dsolve('Dy1=y2*y3','Dy2=-y1*y3','Dy3=-0.51*y1*y2','y1(0)=0,y2(0)=1,y3(0)=1','x')
%% 调用ode45函数求解微分方程df2,自变量为x,范围为[0,4*pi] ; 因变量为y1,y2和y3,初始值: y1(0)=0,y2(0)=y3(0)=1
[x, y] = ode45('df2', [0 4*pi], [0 1 1]); % 这里的y是一个有3列的矩阵哦!
function dy = df2(x,y)
% 注意哦,x是自变量,y是因变量,由y1,y2,y3组成
dy = zeros(3,1); % 初始化用来储存因变量一阶导数的列向量(不能写成行向量哦)
dy(1) = y(2) * y(3);
dy(2) = -y(1) * y(3);
dy(3) = -0.51 * y(1) * y(2);
% 上面四行可以写成一行: dy = [ y(2) * y(3); -y(1) * y(3); -0.51 * y(1) * y(2)]
en
(3)二阶微分方程
% 我们先看这个微分方程有没有解析解
dsolve('(1+x*x)*D2y=2*x*Dy','y(-2)=3,Dy(-2)=4','x')
[x,y]=ode45('df3',[-2,2],[3,4]); % 求出这个微分方程df3的数值解
function dy=df3(x,y)
dy=zeros(2,1); % 初始化用来储存因变量一阶导数的矩阵
dy(1)=y(2);
dy(2)=(2*x)/(1+x*x)*y(2);
en