MATLAB教学_11方程式求根

本文视频地址：https://www.bilibili.com/video/av68228488?p=14

本文学习三个内容：

1. Symbolic approach
2. Numeric root solvers
3. Rescursive functions

目录

Problem Statement

sym()

syms

solve()

10分钟练习

12分钟练习

21分钟练习

diff()

int()

fsolve()

fzero()

optimset

roots()

Numeric Root Finding Methods

迭代函数

---

Problem Statement

求方程的根

1. 解析解。用的以前学的公式。 
2. 图解积分法。用图来看为0的值在哪
3. 数值解析

sym()

创建符号变量，表达式，函数，矩阵。

要创建符号表达式，首先要创建符号变量，然后对它们进行操作。

参考地址 ：https://ww2.mathworks.cn/help/symbolic/sym.html

syms

syms var1 ... varN

创建符号变量var1…varN。也就是把var1...varN，设置为符号变量，可以用到求解方程式里。用空格分隔不同的变量。syms从变量中清除了所有的假设。

syms 和 sym() 功能基本相同，只是表达方式不同；

%%
syms x    %5分钟练习
x+x+x
(x+x+x)/4
y=sym('y');
y+y+y
(y+y+y)/4

solve()

方程组求解器函数。

S = solve(eqn,var) 

为变量var求解方程eqn。如果不指定var，则symvar函数确定要求解的变量。例如，解(x + 1 == 2, x)解出x = 2的方程。

S = solve(eqn,var,Name,Value) 

这个格式里面 'Name' 是属性的设置。有关属性的设置可以参考网址：

https://ww2.mathworks.cn/help/symbolic/solve.html#bue1f92-1

输出一般用 S 来表示：也就是方程组的解。

方程的解，作为符号数组返回。符号数组的大小对应于解的数目。

 

10分钟练习

clear all;  clc;   %10分钟练习
syms x;
y=x*sin(x)-x;
S=solve(y,x)     %方程的解一般用S表示.
Z=solve(x*sin(x)-x,x)   %如果方程没有右端的话，默认是=0
solve(cos(x)^2-sin(x)^2,x)
solve(cos(x)^2+sin(x)^2,x)

12分钟练习

%%
clear all;  clc;   %12分钟练习
syms x y
eq1=x-2*y-5;
eq2=x+y-6;
A=solve(eq1,eq2,x,y)    %放在A里，A相当于一个struct 结构
[X Y]=solve(eq1,eq2,x,y)   %这样可以分别显示 x 和 y 的解。 和函数里  x,y的顺序一样

clear all;  clc;   %12分钟练习
syms a b x
S=solve(a*x^2-b,x)

21分钟练习

%%
clear all;  clc;   %21分钟练习
syms x y a b r
solve((x-a)^2-(y-b)^2==r^2,x)
%%
clear all;  clc;   %21分钟练习
syms a b c d
syms x y z g
A=[a b;c d]*[x y;z g];    %构建求逆矩阵的公式
B=[a b;c d]
inv(B)    %先用函数inv() 算下，是否正确
[X Y Z G]=solve(A==eye(2),x,y,z,g);   %构建公式，
C=[X Y; Z G]   %显示结果，与inv()所算相同

diff()

微分符号表达式或函数

diff(F,var)   %对变量var进行微分。

diff()发现这里有两个用法。下面是在MATLAB教学10里面的用法。是求差分的。

Y = diff(X) 计算沿大小不等于 1 的第一个数组维度的 X 相邻元素之间的差分：

• 如果 X 是长度为 m 的向量，则 Y = diff(X) 返回长度为 m-1 的向量。Y 的元素是 X 相邻元素之间的差分。

  Y = [X(2)-X(1) X(3)-X(2) ... X(m)-X(m-1)]

clear all;  clc;   %23分钟练习
syms x y
z=4*x^5;
zprime=diff(z)
f=exp(x^2)/(x^3-x+3);
fprime=diff(f)
g=(x^2+x*y-1)/(y^3+x+3);
gprime=diff(g,x)      %因为这里有两个未知数，所以指导为x

int()

定积分和不定积分。

F = int(expr) 

计算expr的不定积分。int使用由symvar(expr,1)确定的默认集成变量。如果expr是常量，那么默认的积分变量是x。

F = int(expr,var)

计算expr对符号标量变量var的不定积分。

F = int(expr,a,b) 

计算从a到b的定积分。int使用由symvar(expr,1)确定的默认积分变量。如果expr是常量，那么默认的积分变量是x。
int(expr，[a b])等价于int(expr,a,b)。

symvar(expr,1)

指的是在expr 函数，返回1个离x 最近的这个变量。比较难理解，见下面网址：

https://ww2.mathworks.cn/help/symbolic/symvar.html#mw_561e9790-6995-44e9-9ae1-974d1bf3eb87

subs(s,old,new)

返回s的副本，将 s 中的 old替换为new,然后计算s。

clear all;  clc;   %25分钟练习
syms x
y=x^2*exp(x);
z=int(y)   %求 y 的不定积分
z=z-subs(z,x,0)  %这里 subs(),也就是求了z(0)

这里我们先求出了z,然后发现 z(0) 不等0。而是等于2。那么根据积分，我们知道肯定是有一个常数是2。所以用 subs(z,x,0)来求出这个2。

%%
clear all;  clc;   %28分钟练习
syms x;
y=(x^2-x+1)/(x+3);
z=int(y,[0 10])

fsolve()

解非线性方程组。

非线性系统解算器，解决指定的问题F (x) = 0。是针对于 x 来说， 会返回一个向量或矩阵。

x = fsolve(fun,x0) 

从x0开始，尝试解出有趣的(x) = 0的方程，这是一个由0组成的数组。那么先假设fun(x0)=0,然后在x0周围解这个函数。解出来的值 就是 f(x)=0

%%
clear all;  clc;   %33分钟练习
f2=@(x) (1.2*x+0.3+x*sin(x));
fsolve(f2,0)     %先假设令 f2(x)=0 的 x=0，在这个周围求解这个函数，返回的是真正让这个函数为0的值。

%%
clear all;  clc;   %35分钟练习
f{1}=@(x) (2*x(1)-x(2)-exp(-x(1)));   %错误，这里f()不行，只能用f{}。
f{2}=@(x) (-x(1)+2*x(2)-exp(-x(2)));
fsolve(f,[-5,-5])    %算出来答案有问题，不知道哪里有错。

%%
clear all;  clc;   %35分钟练习
f=@(x) [2.*x(1)-x(2)-exp(-x(1));-x(1)+2.*x(2)-exp(-x(2))];   %参照弹幕大神写的，是把f 设为一个矩阵。
fsolve(f,[-5,-5])

还有一个办法是建立一个函数。然后直接调用这个函数。这个函数要以函数为名单独建立一个文件。

fzero()

 

x = fzero(fun,x0) 尝试求出 fun(x) = 0 的点 x。此解是 fun(x) 变号的位置 - fzero 无法求函数（例如 x^2）的根。

也就是说fun(x) 必须要穿过x 轴，这样fun(x)才会由正变为负，或者由负变为正。

clear all;  clc;   %37分钟练习
f=@(x) x.^2;
fzero(f,0.1)    %l因为x^2没有穿过x轴，所以用fzero()解不出来
fsolve(f,0)

optimset

创建或修改优化 options 结构体。主要是对options进行修改，里面有若干参数，具体可见以下网：

https://ww2.mathworks.cn/help/matlab/ref/optimset.html

roots()

roots() 只能用来算多项式

%%
clear all;  clc;    %43分钟练习
roots([1 -3.5 2.75 2.125 -3.875 1.25])

%%
clear all;  clc;    %44分钟练习
roots([1 -6 -12 81])

 

Numeric Root Finding Methods

第一种方法就是用包夹的区间，越来越小。

第二种是牛顿迭代法；沿着导直到接近0。

什么时候计算结束：用精度 或者是 迭代的次数

 

 

这个使用时函数要有单调性，且是连续的。负责会出错。

大概思路就是取这段区间两头的值，那么 相乘肯定是负值，说明区间取值肯定把 f(x)= 0 的 x 值包进去。再取小一点的值 ，一直到区间赵接近 f(x)=0

 

牛顿迭代法，先假设一个点，然后求导与x 轴相交，作为新的点，再求导，直到与x 轴相交的 f(x) 的点找到。

这是得到 x 的值。 先假设为 x0,那么根据斜率， slope=f(x0)/a =f'(x0),那么 a=f(x0)/f'(x0), 这样长度就知道了。

所以 x1=x0-a=x0-f(x0)/f'(x0)

迭代函数

function output=fact(n)    %自建的函数，可以直接调用
if n==1
    output=1;      %如果 n=1,则直接输出为1
else
    output=n*fact(n-1);   %如果n >1,则一直迭代下去
end
end