代数方程和微分方程求解

非线性方程组的一般形式为
F(x)=0
其中
f1 ( x1 , x 2 ,, x n ) f 2 ( x1 , x 2 ,, x n ) F ( x) f ( x , x ,, x ) n n 1 2
非线性方程组求解的matlab函数为 x=fsolve(fun,x0); 使用方法完全类似于函数fzero，但对初始条 件的选择要困难的多。 例4.9：求解非线性方程组
例4.7：求方程 的最小正根。
kx sin(kx) 10
k 1
5
首先寻找初始值。画出方程
f ( x) kx sin(kx) 10
k 1
5
在[0,2] 上的图形 >> x=0:pi/50:2*pi; >> y=x.*sin(x)+2*x.*sin(2*x)+3*x.*sin(3*x)... +4*x.*sin(4*x)+5*x.*sin(5*x)-10; >> plot(x,y,x,0)

solve函数还有另一种调用形式 x=solve(eq1,eq2,…,eqm,var1,…,varm); 其中 eq 是字符串，定义方程。 var也是字符串， 指明所求解的变量。
例4.1 求解线性方程组
ax y z 1 x ay z 2 x y az 3
参考程序： n=1000; A=diag(2*ones(1,n))+diag(ones(1,n1),1)+diag(ones(1,n-1),-1); A(n,1)=1;b=ones(n,1); tic x=A\b; toc A1=sparse(A); tic x=A1\b; toc 输出结果： Elapsed time is 0.144504 seconds. Elapsed time is 0.009240 seconds.
例4.5：求方程 sin(x)+cos(x)=1 的解。初始值依次选取 x0=1,2,3,…,10 。
参考程序： r=zeros(1,10); for x0=1:10 r(x0)=fzero(@(x)(sin(x)+cos(x)-1),x0); end r 运行结果： r = 1.5708 1.5708 1.5708 1.5708 6.2832 6.2832 6.2832 7.8540 7.8540 7.8540 运行结果显示：在一般情况下，迭代收敛到与初始值接 近的根。
输出结果：
由上图知，方程组共有20个不同的解。再由方 程组关于x的对称性知只要求出右边的10个。利 用局部放大图（如下图）可得到10个解的近似 值
数值程序： function r=d2solve h=ezplot(@(x,y)(x.^2.*cos(pi*x)+y.^2.*sin (pi*y)-pi/2),[-3,3,-3,3]) str1=get(get(gca,'Title'),'String') set(h,'Color',[1 0 0],'LineStyle',':') hold on h1=ezplot(@(x,y)(x.^2+y.^2+2*sin(2*x.^2.* y).^2-4),[-3,3,-3,3]) str2=get(get(gca,'Title'),'String') title([str1,',',str2])
Y=polyval(p,x)
的值
% 输出多项式 p 在向量 x
多项式求根的函数为 r=roots(p) 求得多项式p的所有根。 例4.3：求多项式
x 6 20x 5 138x 4 328x 3 223x 2 1692x 1260
的根并在多项式图形中表示。
参考程序：
q =[1 -20 138 -328 -223 1692 1260]; r=roots(q); x=-2.2:0.05:8; y=polyval(q,x); y1=polyval(q,r); plot(x,y,r,y1,'p') xlim([-2.2,8]) legend('polynomial','roots')
2 2 x cos(x) y sin(y ) 2 2 2 2 2 x y 2 sin( 2 x y ) 4
的全部实根。
解：首先将两条联立的曲线画出，以确定根（交 点）的个数和位置。 参考程序： h=ezplot(@(x,y)(x.^2.*cos(pi*x)+y.^2.*... sin(pi*y)-pi/2),[-3,3,-3,3]) str1=get(get(gca,'Title'),'String') set(h,'Color',[1 0 0],'LineStyle',':') hold on h1=ezplot(@(x,y)(x.^2+y.^2+2*sin(2*x.^... 2.*y).^2-4),[-3,3,-3,3]) str2=get(get(gca,'Title'),'String') title([str1,',',str2]) grid on
程序运行结果 >> r=d2solve r = 1.7251 1.7038 1.6679 1.6430 0.8980 0.8899 -0.6279 -0.4121 0.4197 0.7085 1.7007 -1.5239
Biblioteka Baidu
利用对上述图形的局部放大图（下图）估计得到初始 值3.8
程序和结果： >> format long >> r=fzero(@(x)(x.*sin(x)+2*x.*sin(2*x)+3*x .*sin(3*x)+4*x.*sin(4*x)+5*x.*sin(5*x)10),3.8) r = 3.908485977246559
例4.8：求含参数代数方程
f ( x, a) a 3x e
2 0.2 x
e
0.02 x
sin ax
4
的根与参数a的关系。
注：含参数的方程的求解对后续计算非常重要。 如上面例子求参数，使得根取极大等。
function rootfinding x=0:0.01:5;l=1;root_1=zeros(1,length(x)); for a=x; root_1(l)=fzero(@(x)fun(x,a),a); l=l+1; end plot(x,root_1,'-') xlabel('parameters') ylabel('roots') end function f=fun(x,a) f=a-3*x.^2.*exp(-0.2*x)+exp(-0.02*x).*sin(a*x).^4; end
代数方程组的符号求解函数为 solve 函数，函数 调用格式为 x=solve(fun1,fun2,…,funm,var1,…,varm) 其中fun1, fun2,…为符号表达式。 var1,var2,….为符号变量，指明所求解的变量。 如果不指明所求解的变量，则系统根据变量在英 文字母中与x的接近程度选择。
参考程序 syms x y z a [u,v,w]=solve(a*x+y+z-1,x+a*y+z-... 2,x+y+a*z-3,x,y,z) 输出结果： u = (a - 4)/(a^2 + a - 2) v = 2/(a + 2) w = (3*a)/(a^2 + a - 2)
也可以利用下面的语句求解 >> s=solve('a*x+y+z-1','x+a*y+z-2','x+y+a*z3','x','y','z') s = x: [1x1 sym] y: [1x1 sym] z: [1x1 sym] 其中，s是结构形数据。可以用下面的方法显示求得的解 >> s.x ans = (a - 4)/(a^2 + a - 2) 方程组的符号求解得到的是解析式。其优点是可以得到含 未知参数的解，有利于进一步的分析。缺点是计算效率低， 且复杂的问题一般得不到解析解。
例4.6：求解方程
x e dt 1 0
t 2 0
x
function r=rootexample(x0) r=fzero(@rootfun,x0); end function f=rootfun(x) f=x+quad(@(t)exp(-t.^2),0,x)-1; end 运行结果： >> r=rootexample(0.8) r = 0.5219
求解非线性方程 f(x)=0 一般利用以下形式的迭代方法： xk+1=(xk) k=1,2,… 如著名的Newton迭代法

xk 1
f ( xk ) f ( xk )
就是众多迭代方法中的一种。
迭代是一种逐步近似的过程，在迭代求根中，需 要给定初始近似 来启动迭代过程。初始近似 的 选取非常重要，这是因为 1.许多方法只有当 接近所求的解时，迭代才收敛。 2.当方程有多个解时， 需要充分靠近所求的解， 才能保证迭代收敛到这个解。 Matlab 中非线性方程求根的函数是 fzero ，调用格 式为 p=fzero(fun,x0); 其中x0是初始值。fun可以是字符串以指示函数名， 也可以是以@开头的函数句柄如 ‘sin’，@sin，@(x)sin(x) 都可以，应用最灵活的是后一种。
x=[1.72,1.7,1.68,1.68,1.66,1.64,0.36,0.65,0.9,0.88]; y=[-0.63,-0.42,-0.19,0.16,0.42,0.7,-1.85,-1.38,1.7,-1.53]; X=[x;y];m=length(x); r=zeros(2,m); for k=1:m r(:,k)=fsolve(@solfun,X(:,k)); end plot(r(1,:),r(2,:),'ko') end function f=solfun(x) f=[x(1).^2.*cos(pi*x(1))+x(2).^2.*sin(pi*x(2))pi/2 x(1).^2+x(2).^2+2*sin(2*x(1).^2.*x(2)).^2-4]; end

例4.2：求微分方程
d2y dy 2 2y 0 2 dt dt 的通解和在初始条件y(0)=1,y’(0)=2下的特解
>> y=dsolve('D2y-2*Dy+2*y=0') y = C24*exp(t)*cos(t) + C25*exp(t)*sin(t) >>y=dsolve('D2y-2*Dy+2*y=0',... 'y(0)=1','Dy(0)=2') y = exp(t)*cos(t) + exp(t)*sin(t)
-

线性方程组 Ax=b 可以利用矩阵除法直接得到。但当系数矩阵为稀疏 矩阵时，利用稀疏矩阵函数可以得到更高的计算效 率。 稀疏矩阵利用函数 A1=sparse(A); 定义。
例4.4：求解n阶线性方程组
0 x1 1 2 1 1 2 1 x 2 1 1 2 x3 1 1 0 x 1 0 1 2 n 1 对n=1000，分别利用正常方法和稀疏矩阵求 解，并比较计算时间。
常微分方程（组）的符号运算函数为dsolve，调 用格式为 变量=dsolve(eq1,eq2,…); 其中equ1, equ2等是字符串，描述微分方程。方 程形式与solve语句中类似，其中的导数项利用D 表示， Dy,D2y 分别描述 y 的一阶和二阶导数。初 值和边值条件也利用方程描述。
多项式求根
多项式 在matlab语言中，多项式利用行向量表示。如向量 u=[1 2 0 -5 4] 用作多项式函数时，与多项式

p x 4 2 x 3 5x 4
对应。因此，多项式运算对应向量的相关运算，如 多项式的加法对应向量的加法等。
多项式运算的几个常用函数： P=conv(p1,p2); [d,r]=deconv(p1,p2); Dp=polyder(p); Ip=polyint(p) %多项式乘法 %多项式除法 %多项式的导数 %多项式的积分（原函数）