实验8 非线形方程和常微分方程的解法Word下载.docx
《实验8 非线形方程和常微分方程的解法Word下载.docx》由会员分享，可在线阅读，更多相关《实验8 非线形方程和常微分方程的解法Word下载.docx（9页珍藏版）》请在冰豆网上搜索。
得出结果：
ans=
0.50183.14076.28329.4248
【例】求解方程组x－0.7sinx－0.2cosy=0
y－0.7cosx+0.2siny=0
先编制函数文件fu.m：
functiony=fu（x）
y
（1）=x
（1）-0.7*sin（x
（1））-0.2*cos（x
（2））;
y
（2）=x
（2）-0.7*cos（x
（1））+0.2*sin（x
（2））;
y=[y
（1）,y
（2）];
在命令窗口调用fu运算：
x1=fsolve（'
fu'
[0.5,0.5]）
0.52650.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
注意：
1.方程解的估计值x0可用fplot作图看出；
2.用function建立函数文件fn，求解调用时fn两边要加单引号，而用inline时fn两边不要加单引号；
3.这两种方法也可解线性方程组.
2.代数方程的符号解
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.
【例】
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'
）
计算结果为：
[1/2/a*（-b+（b^2-4*a*c）^（1/2））]
[1/2/a*（-b-（b^2-4*a*c）^（1/2））]
-（a*x^2+c）/x
x=
4/3
y=
-1/3
a=
[2]
[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.28中的方程，结果不对，所以要验根，多用几种方法相互验证.用符号法解方程3x2－ex=0，解的表达式不易懂，怎么办？
>
x=solve（'
3*x^2－exp（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）
[.912]
[3.72]
[-.460]
3.常微分方程数值解法（ordinarydifferentialequation）
[T,Y]=solver（odefun,tspan,y0）%在区间tspan=[t0,tf]上，用初始条件y0求解显式微分方程y'
=f（t,y）.
1.solver为命令ode45，ode23，…之一.
2.odefun为显式常微分方程y'
3.tspan积分区间（即求解区间）的向量tspan=[t0,tf].要获得问题在其他指定时间点t0,t1,t2,…上的解，则令tspan=[t0,t1,t2,…,tf]（要求是单调的）.
4.y0包含初始条件的向量.
求解具体ODE的基本过程如下所示.
1.用微分方程与初始条件描述所给问题.
F（y,y'
y"
…,y（n）,t）=0
y（0）=y0，y'
（0）=y1，…，y（n－1）（0）=yn－1
而y=[y;
y
（1）;
y
（2）;
…；
y（m－1）]，n与m可以不等.
2．运用数学中的变量替换：
yn=y（n－1），yn－1=y（n－2），…，y2=y'
，y1=y，把高阶（2阶及以上）的方程（组）写成一阶微分方程组：
，
3．根据1与2的结果，编写能计算导数的M函数文件odefile.
4．将文件odefile与初始条件传递给求解器Solver中的一个，运行后就可得到解列向量y
因为没有一种算法可以有效地解决所有的ODE问题，为此，MATLAB提供了多种求解器Solver，对于不同的ODE问题，采用不同的Solver，具体如表1.11所示.
表1.11不同求解器Solver的特点
求解器Solver
ODE类型
特点
说明
ode45
非刚性
一步算法，4，5阶Runge-Kutta方程，累计截断误差达（△x）3
大部分场合的首选算法
ode23
一步算法，2，3阶Runge-Kutta方程，累计截断误差达（△x）3
使用于精度较低的情形
ode113
多步法，Adams算法，高低精度均可到10－3～10－6
计算时间比ode45短
ode23t
适度刚性
采用梯形算法
适度刚性情形
ode15s
刚性
多步法，Gear’s反向数值微分，精度中等
若ode45失效时，可尝试使用
ode23s
一步法，2阶Rosebrock算法，低精度
当精度较低时，计算时间比ode15s短
【例】求解微分方程y'
=－2y+2x2+2x，0≤x≤0.5，y（0）=1.
-2*y+2*x^2+2*x'
x'
y'
[x,y]=ode23（fun,[0,0.5],1）;
x'
Columns1through7
00.04000.09000.14000.19000.24000.2900
Columns8through12
0.34000.39000.44000.49000.5000
y'
1.00000.92470.84340.77540.71990.67640.6440
0.62220.61050.60840.61540.6179
图形结果如图1.14所示.
图1.14例1.76图形结果图1.15例1.77图形结果
【例】求解描述振荡器的经典的VerderPol微分方程
，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：
functionfy=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）
图形结果如图1.15所示.
4.常微分方程的符号解
r=dsolve（‘eq1,eq2,…’,‘cond1,cond2,…’,‘v’）.
1.对给定的常微分方程（组）eq1，eq2，…中指定的符号自变量v，与给定的边界条件和初始条件cond1，cond2，…求符号解（即解析解）r；
2.若没有指定变量v，则缺省变量为t；
在微分方程（组）的表达式eq中，Dy=dy/dx，D2y=d2y/dx2；
3.初始和边界条件由字符串表示：
y（a）=b，Dy（c）=d，D2y（e）=f，等等，分别表示
；
4.若边界条件少于方程（组）的阶数，则返回的结果r中会出现任意常数C1，C2.
【例】
D1=dsolve（'
D2y=Dy+exp（x）'
D2=dsolve（'
（Dy）^2+y^2=1'
s'
D3=dsolve（'
Dy=a*y'
y（0）=b'
）%带一个初始条件
D4=dsolve（'
D2y=-a^2*y'
y（0）=1'
Dy（pi/a）=0'
）%带两个初始条件
[x,y]=dsolve（'
Dx=y'
Dy=-x'
）%求解线性微分方程组
D1=
-exp（x）*t+C1+C2*exp（t）
D2=
[-1]
[1]
[sin（s-C1）]
[-sin（s-C1）]
D3=
b*exp（a*t）
D4=
cos（a*t）
cos（t）*C1+sin（t）*C2
-sin（t）*C1+cos（t）*C2
三、练习和思考
①求解方程3x2－ex=0
②求解方程3x2－lnx=10
③求解方程5x4－sin2x=0
④求解微分方程y"
－y'
+2y=5，y（0）=1，y'
（0）=2.
⑤求解微分方程的特解，并作出解函数曲线图.
y"
－2（1－y2）y'
+y=0，0≤x≤30，y（0）=1，y'
（0）=0
升级会员 