实验8 非线形方程和常微分方程的解法.docx
- 文档编号:9129752
- 上传时间:2023-02-03
- 格式:DOCX
- 页数:9
- 大小:66.34KB
实验8 非线形方程和常微分方程的解法.docx
《实验8 非线形方程和常微分方程的解法.docx》由会员分享,可在线阅读,更多相关《实验8 非线形方程和常微分方程的解法.docx(9页珍藏版)》请在冰豆网上搜索。
实验8非线形方程和常微分方程的解法
实验8非线性方程和常微分方程的解法
一、实验目的
求解非线性方程和常微分方程.
二、实验内容与要求
1.非线性方程的数值解
(1)最小二乘法
格式:
fsolve(‘fun’,x0)%求方程fun=0在估计值x0附近的近似解.
【例】求方程x-e-x=0的解.
>>fc=inline('x-exp(-x)');
>>x1=fsolve(fc,0)
x1=
0.5671
f(x)=5x2sinx-e-x图像
问题:
求解方程5x2sinx-e-x=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],1e-6)
得出结果:
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])
x1=
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')
>>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.28中的方程,结果不对,所以要验根,多用几种方法相互验证.用符号法解方程3x2-ex=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.常微分方程数值解法(ordinarydifferentialequation)
格式:
[T,Y]=solver(odefun,tspan,y0)%在区间tspan=[t0,tf]上,用初始条件y0求解显式微分方程y'=f(t,y).
说明:
1.solver为命令ode45,ode23,…之一.
2.odefun为显式常微分方程y'=f(t,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.
>>fun=inline('-2*y+2*x^2+2*x','x','y');
>>[x,y]=ode23(fun,[0,0.5],1);
>>x'
ans=
Columns1through7
00.04000.09000.14000.19000.24000.2900
Columns8through12
0.34000.39000.44000.49000.5000
>>y'
ans=
Columns1through7
1.00000.92470.84340.77540.71990.67640.6440
Columns8through12
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)
x=
cos(t)*C1+sin(t)*C2
y=
-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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 实验8 非线形方程和常微分方程的解法 实验 线形 方程 微分方程 解法