非线性方程和常微分方程的解法

非线性方程和常微分方程的解法

实验8 非线性方程和常微分方程的解法

一、实验目的

学会用MATLAB软件求解非线性方程和常微分方程。

二、实验内容与要求

1. 非线性方程的整值解

（1）最小二乘法

格式：fsolve（’fun’,x0）%求方程fun=0在估计值x附近的近似解。

[例1.72] 求方程x e 0的解。

>>fc=inline(‘x-exp(-x)’);

>>xl=fsolve(fc,0)

xl=

0.5671

问题1.28：求解方程5xsinx-e 0,观察知有多解，如何求之？

先用命令fplot（’5*x^2*sin(x)-exp(-x),0]’,[0,10]）作图1.13，注意5*x^2*sin(x)-exp(-x),0]中的“[ ，0]”是作y=0直线，即x轴。方程在[0，10]区间从图中可看出有4个解，分别在0，3，6，9附近，所以用命令：

>>fun=inline(’5*x^2*sin(x)-exp(-x)’);

>>fsolve(fun,[0,3,6,9],le-6)

得出结果：

ans=

0.5918 3.1407 6.2832 9.4248

【例 1.73】求解方程组x-0.7sinx-0.2cosy=0 2-x-x

y 0.7cosx 0.2siny 0

先编制函数文件 fu.m:

function y=fu (x)

y (1)=x (1) - 0.7*sin ( x (1) ) - 0.2*cos( x (2) ) ;

y (2)=x (1) - 0.7*cos (x (1) ) + 0.2*sin (x (2) );

y =[ y(1), y(2) ];

在命令窗口调用fu运算：

>>x1 =fsolve( ‘fu’, [0.5,0.5 ])

x1 =

0.5265 0.5079

(2) 零点法

格式：fzero('fun',x0) %求函数fun在x0附近的零点。

说明：估计值x0若为标量时，则在x0附近查找零点，x0=[x1,x2]向量时，则首先要求函数值 fun(x1)fun(x2) 0，然后将严格在[x1,x2]区间内零点，若找不到，系统将给出提示。

非线性方程和常微分方程的解法

【例 1.74】 求函数 f(x) sinx2/x xex 4 的零点。

>> fn =inline('sin (x^2) / x+x* exp (x) - 4' );

>>x=fzero ( fn, [1,2] ) %这里的fn不要加单引号

x =

1.0748

注意：方程解的估计值 可用fplot作图看出；用function建立函数文件fn，求解调用时

fn两边要加单引号，而用inline时fn两边不要加单引号；这两种方法也可解线性方程组。

⒉ 代数方程的符号解

格式：g solve(eq) %求解方程eq 0，输入参量eq可是符号表达式或字符表达式。

g solve(eq,var) %对eq中指定的变量var求解方程eq(var) 0。

g solve(eq1,eq2, ,eqn) %求解方程组eq1 0,eq2 0, eqn 0。

g solve(eq1,eq2, eqn,var1,var2, ,varn) %对方程组eq1,eq2, eqn中指定的n

个变量加var1,var2, varn求解

【例1.75】

>>solve ( 'a*x^2 + b*x + c')

>>solve( 'a*x^2 + b*x + c', 'b')

>>[x,y] =solve ( 'x+y = 1', 'x -11*y= 5')

>>[a,u,v] =solve ('a*u^2 + v^2', 'u- v =1','a^2 -5*a +6')

计算结果为：

ans =

[ 1/2/a*(-b +(b^2-4*a*c)^(1/2) ) ]

[ 1/2/a*(-b -(b^2 -4*a*c)^(1/2) )]

ans =

-(a*x^2+c) /x

x =

4/3

y =

-1/3

a =

[2]

[2]

[3]

[3]

u =

[ 1/3+1/3*i*2^(1/2) ]

[ 1/3-1/3*i*2^(1/2) ]

[ 1/4+1/4*i*3^(1/2) ]

非线性方程和常微分方程的解法

[ 1/4-1/4*i*3^(1/2) ]

v =

[ -2/3+1/3*i*2^(1/2)]

[ -2/3-1/3*i*2^(1/2)]

[ -3/4+1/4*i*3^(1/2)]

[ -3/4-1/4*i*3^(1/2)]

注意：对于单个的方程或方程组，若不存在符号解，则返回方程（组）的数值解。

问题1.29：用符号法求解问题1.28中的方程，结果不对，所以要验根，多用几种方法相互验证，用符号法解方程3x e 0，解的表达式不易懂，怎么办？

x =solve('3*x^2-exp(x)')

x =

[ -2*lambertw (-1/6*3^(1/2) ) ]

[ -2*lambertw (-1,-1/6*3^(1/2) ) ]

[ -2*lambertw (1/6*3^(1/2) ) ]

再用命令：

>>vpa (x,3)

ans =

[ .912]

[ 3.72]

[ -.460]

3.常微分方程数值解法

格式：[T,Y] solver(odefun,tspan,y0) %区间tspan [t0,tf]上，用初始条件y0求解

显示微分方程y' f(t,y)

说明：solver为命令ode45,ode23,ode113,ode15s,ode23s,ode23t,ode23tb,之一。

odefun为显示常微分方程y' f(t,y)。

tspan 积分区间 （即求解区间） 的向量tspan [t0,tf]。要获得问题在其他指的时 间点t0,t1,t2 上的解，则令tspan [t0,t1,t2 tf] （要求是单调的）。 2x

y0包含初始条件的向量。

求解具体ODE的基本过程如下所示。

1 根据问题所属学科中的规律、定律、公式，用微分方程与初始条件进行描述。 ○(n)F（y, y’, y’’, , y,t）=0 (n-1)y(0)=y0, y’(0)=y1, ,y(0)=yn-1

而y=[y; y(1);y(2); ;y(m-1)],n 与m可以不等。 (n-1)(n-2)2 运用数学中的变量替换：yn=y,yn-1=y, ,y2=y’,y1=y,把高阶（大于2阶）的○

方程（组）写成一阶微分方程组：

非线性方程和常微分方程的解法

y1 0 y0 y1' f1 t,y ' ft,yy0 2 y1 y2 2 ,y . . y' . .0 . . . . n-1 yn' ft,yy0y n n

3 根据○1与○2的结果，编写能计算导数的M函数文件odefile. ○

4 将文件odefile与初始条件传递给求解器Solver中的一个，运行后就可得到ODE的、○

在指定时间区间上的解列向量y（其中包含y及不同阶的导数）。

因为没有一种算法可以有效地解决所有的ODE问题，为此，MATLAB提供了多种求解器Solver，对于不同的ODE问题，采用不同的Solver,具体如表1.11所示。

【例1.76】 求解微分方程y’=－2y+2x+2x, 0≦x≦0.5, y(0)=1.

>> fun=inline('-2*y+2*x^2+2*x','x','y');

>> [x,y]=ode23(fun,[0,0.5],1);

>> x'

ans =

Columns 1 through 6

0 0.0400 0.0900 0.1400 0.1900 0.2400

Columns 7 through 12

0.2900 0.3400 0.3900 0.4400 0.4900 0.5000

>> y'

ans =

Columns 1 through 7

1.0000 0.9247 0.8434 0.7754 0.7199 0.6764 0.6440

Columns 8 through 12

0.6222 0.6105 0.6084 0.6154 0.6179

【例１．７７】求解描述振荡器的经典的Ver der Pol 微分方程

d2y2dy'2

dt2 (1 y)dt y 0,y(0) 1,y(0) 0.

分析：令

非线性方程和常微分方程的解法

x1 y,x2 dy/dt, 7,则：

dx1/dt x2

dx2/dt 7(1 x1^2)x2 x1

编写函数文件vdp,m:

function fy=vdp(t,x)

fy=[x(2);7*(1-x(1)^2)*x(2)-x(1)];

再在命令窗口中执行：

>> Y0=[1;0];

>> [t,x]=ode45('vdp',[0,40],Y0);

>> y=x(:,1);dy=x(:,2);

>> plot(t,y,t,dy)

图形结果如图１.１５所示

.

图１.１５ 例１.７７图形结果

４．常微分方程 …… 此处隐藏：2853字，全部文档内容请下载后查看。喜欢就下载吧 ……