例子.docx
- 文档编号:17823813
- 上传时间:2023-08-04
- 格式:DOCX
- 页数:26
- 大小:305.51KB
例子.docx
《例子.docx》由会员分享,可在线阅读,更多相关《例子.docx(26页珍藏版)》请在冰点文库上搜索。
例子
下载MATLAB的工程应用:
实验5符号方程的求解
MATLAB7.0中的符号计算可以求解线性方程(组)、代数方程的符号解、非线性符号方程(组)、常微分方程(组),求解这些方程(组)是通过调用solve函数实现的,如求解代数方程的符号解调用solve函数的格式是solve('eq')、solve('eq','v')、[x1,x2,…xn]=solve('eq1','eq2',…'eqn')等,求解非线性符号方程是调用优化工具箱的fsolve函数,调用格式有fsolve(f,x0)、fsolve(f,x0,options)、[x,fv]=fsolve(f,x0,options,p1,p2…)等,而解常微分方程(组)则是调用dsolve函数,调用的格式有[x1,x2,…]=dsolve('eq1,eq2,…','cond1,cond2…','v')。
现将各函数的调用格式列于下表(表5—1),在各个实例中说明各种格式的用法。
表5—1符号方程求解的solve函数调用格式
调用格式
说明
solve('eq')
对系统默认的符号变量求方程eq=0的根。
solve('eq','v')
对指定变量v求解方程eq(v)=0的根。
[x1,x2,…xn]=solve('eq1','eq2',…'eqn')
对系统默认的一组符号变量求方程组eqi=0(i=1,2,…n)的根。
[v1,v2,…vn]=solve('eq1','eq2',…'eqn','v1','v2',…'vn')
对指定的一组符号变量v1,v2,…vn求方程组eqi=0(i=1,2,…n)的根。
linsolve(A,B)
求符号线性方程(组)AX=B的解。
相当于X=sym(A)\sym(B)
fsolve(f,x0)
从x0开始搜索f=0的解。
fsolve(f,x0,options)
根据指定的优化参数options从x0开始搜索f=0的解。
fsolve(f,x0,options,p1,p2…)
优化参数option不是默认时,在p1,p2…条件下求f=0解。
优化参数option可取的值有0(默认)和1
[x,fv]=fsolve(f,x0,options,p1,p2…)
优化参数option为默认时,在p1,p2…条件下求f=0解,并输出根和目标函数值。
[x,fv,ex]=fsolve(f,x0,options,p1,p2…)
优化参数option为默认时,在p1,p2…条件下求f=0解,并输出根和目标函数值,并通过exitflag返回函数的退出状态。
[x,fv,ex,out]=fsolve(f,x0,options,p1,p2…)
优化参数option为默认时,在p1,p2…条件下求f=0解,并给出优化信息。
[x,fv,ex,out,jac]=fsolve(f,x0,options,p1,p2…)
优化参数option为默认时,在p1,p2…条件下求f=0解,输出值为x处的jacobian函数。
[x1,x2,…]=dsolve('eq1,eq2,…','cond1,cond2…','v')
在初始条件为cond1,cond2…时求微分方程组eq1,eq2,…对指定变量v的特解。
[x1,x2,…]=dsolve('eq1','eq2',…,'cond1','cond2'…,'v')
同[x1,x2,…]=dsolve('eq1,eq2,…','cond1,cond2…','v')
一、代数方程的符号解
MATLAB7.0中求代数方程的符号解是通过调用solve函数实现的。
用solve函数求解一个代数方程时的调用格式一般是:
solve('代数方程','未知变量')或x=solve('代数方程','未知变量')
当未知变量为系统默认变量时,未知变量的输入可以省略。
当求解由n个代数方程组成的方程组时调用的格式是:
[未知变量组]=solve('代数方程组','未知变量组')
未知变量组中的各变量之间,代数方程组的各方程之间用逗号分隔,如果各未知变量是由系统默认的,则未知变量组的输入可以省略。
实例1、求解高次符号方程
和方程
对y的解。
>>symsxyzab%定义符号变量
>>solve(x^4-3*a*x^2+4*b)%求解高次方程
ans=
1/2*(6*a+2*(9*a^2-16*b)^(1/2))^(1/2)
-1/2*(6*a+2*(9*a^2-16*b)^(1/2))^(1/2)
1/2*(6*a-2*(9*a^2-16*b)^(1/2))^(1/2)
-1/2*(6*a-2*(9*a^2-16*b)^(1/2))^(1/2)
>>solve(x^3+2*a*x*y-3*b*y^2,y)%对指定变量求解方程
ans=
1/6/b*(2*a+2*(a^2+3*b*x)^(1/2))*x
1/6/b*(2*a-2*(a^2+3*b*x)^(1/2))*x
实例2、求解多元高次方程组
>>[x,y]=solve('x^3+2*x*y-3*y^2-2','x^3-3*x*y+y^2+5')%求解多元高次方程组
x=
1.8061893129091900210106914427639+1.1685995398225344682988775209345*i
.51233671712308192620449202726936+1.0694475803263816285960240820218*i
-1.2247760300322719472151834700333+.35066213508454219362158900429401*i
-1.2247760300322719472151834700333-.35066213508454219362158900429401*i
.51233671712308192620449202726936-1.0694475803263816285960240820218*i
1.8061893129091900210106914427639-1.1685995398225344682988775209345*i
y=
1.8086294126483514370835126464657+1.9432962587476317909683476452237*i
.17307087932198664953847299268063-.78620181218420502898925154555661*i
-.61451279197033808662198563914677-.89207785198625780793629825881329*i
-.61451279197033808662198563914677+.89207785198625780793629825881329*i
.17307087932198664953847299268063+.78620181218420502898925154555661*i
1.8086294126483514370835126464657-1.9432962587476317909683476452237*i
实例3、求解方程组
的解。
>>[x,y,z]=solve('x-2*y-4','x^2-2*x*y+y-z','x^2-y*z+z')
x=
29/5-1/5*721^(1/2)
29/5+1/5*721^(1/2)
y=
9/10-1/10*721^(1/2)
9/10+1/10*721^(1/2)
z=
241/10-9/10*721^(1/2)
241/10+9/10*721^(1/2)
实例4、求解超越方程
的解。
>>solve('x*2^x-1')%求解超越方程
ans=
1/log
(2)*lambertw(log
(2))
注:
lambertw是一个函数,lambertw(x)表示方程w*exp(w)=x的解w。
其数值可以在命令窗口输入该函数得到。
>>lambertw(log
(2))
ans=
1.4444
二、符号线性方程(组)的求解
符号线性方程(组)的求解与数值线性方程(组)的求解方法相同,采用矩阵左除或函数linsolve,格式为:
X=A\B或X=sym(A)\sym(B)或X=linsolve(A,B)。
其中A为线性方程组的系数矩阵,B为方程右侧的常数列矩阵。
实例5、求符号线性方程组
的符号解。
>>A=sym('[123;-192;203]');%定义符号矩阵A
>>B=[a;b;1];%定义符号矩阵B
>>x=A\B%求解方程
x=
6/13*b+23/13-27/13*a
3/13*b+5/13-7/13*a
-4/13*b-11/13+18/13*a
三、非线性符号方程的求解
非线性符号方程(组)F(X)=0中X是一个向量,求解显示的结果也是一个向量。
它不仅可以用调用solve函数求解,也可以调用函数fsolve求解,而函数fsolve不是MATLAB符号工具箱的函数,它位于优化工具箱内。
实例6、求解非线性符号方程组
,用solve函数和fsolve函数起始点为x0=[0;0]各自求解。
(1)solve函数求解
>>symsx1x2%定义符号变量
>>[x1,x2]=solve('x1-3*x2=sin(x1)','2*x1+x2=cos(x2)','x1','x2')%求解方程组
x1=
.49662797440907460178544085171994
x2=
.67214622395756734146654770697884e-2
(2)fsolve函数求解
先在文件编辑窗口编写如下M文件,并存于系统的work目录下。
functionF=myfun(x)
F=[x
(1)-3*x
(2)-sin(x
(1));2*x
(1)+x
(2)-cos(x
(2))];
然后在命令窗口求解:
>>x0=[0;0];%设定求解初值
>>options=optimset('Display','iter');%设定优化条件
>>[x,fv]=fsolve(@myfun,x0,options)%优化求解
%MATLAB显示的优化过程
NormofFirst-orderTrust-region
IterationFunc-countf(x)stepoptimalityradius
03121
160.0004233080.50.06171
295.17424e-0100.007514334.55e-0051.25
3129.99174e-0221.15212e-0059.46e-0111.25
Optimizationterminated:
first-orderoptimalityislessthanoptions.TolFun.
x=
0.4966
0.0067
fv=
1.0e-010*
0.3161
0.0018
四、常微分方程的符号解
含有自变量、未知函数和未知函数导数(或微分)的等式叫微分方程。
描述自变量与函数关系的等式叫微分方程的初始条件。
适合微分方程的函数叫微分方程的解。
没有初始条件而求得的解叫微分方程的通解,通解中会包含有与方程阶数相同个数的积分常数C1、C2等;有初始条件且满足初始条件的解叫微分方程的特解,特解一般不含有积分常数。
在MATLAB中,用dsolve函数求解微分方程或微分方程组,dsolve函数参数的输入共有三部分,微分方程、初始条件和自变量。
格式是:
dsolve('微分方程','初始条件','自变量')
微分方程部分的输入与MATLAB符号表达式的输入基本相同,微分或导数的输入是用Dy、D2y、D3y、…来表示y的一阶导数
或
、二阶导数
或
、三阶导数
或
、…。
如果自变量是系统默认的,则自变量输入部分可省略。
dsolve函数的输出部分是该方程(组)的解列表,如果dsolve函数找不到解析解,则系统显示一则错误信息。
实例7、求解微分方程组
在无初始条件和有初始条件
下的解。
(1)无初始条件求解
>>[x,y]=dsolve('D2x+Dy+3*x=cos(2*t)','D2y-4*Dx+3*y=sin(2*t)','t')
x=
1/5*cos(2*t)-1/2*C1*cos(t)+1/2*C2*sin(t)+1/2*C3*cos(3*t)-1/2*C4*sin(3*t)
y=
3/5*sin(2*t)+C1*sin(t)+C2*cos(t)+C3*sin(3*t)+C4*cos(3*t)
(2)有初始条件求解
>>[x,y]=dsolve('D2x+Dy+3*x=cos(2*t)','D2y-4*Dx+3*y=sin(2*t)','Dx(0)=1/5','x(0)=0','Dy(0)=6/5','y(0)=0','t')
x=
1/5*cos(2*t)-3/20*cos(t)+1/20*sin(t)-1/20*cos(3*t)+1/20*sin(3*t)
y=
3/5*sin(2*t)+3/10*sin(t)+1/10*cos(t)-1/10*sin(3*t)-1/10*cos(3*t)
实验十二常微分方程和常微分方程组的求解
一、实验目的:
熟悉Matlab软件中关于求解常微分方程和常微分方程组的各种命令,掌握利用Matlab软件进行常微分方程和常微分方程组的求解。
二、相关知识
在MATLAB中,由函数dsolve()解决常微分方程(组)的求解问题,其具体格式如下:
X=dsolve(‘eqn1’,’eqn2’,…)
函数dsolve用来解符号常微分方程、方程组,如果没有初始条件,则求出通解,如果有初始条件,则求出特解。
例1:
求解常微分方程
的MATLAB程序为:
dsolve('Dy=1/(x+y)','x'),注意,系统缺省的自变量为t,因此这里要把自变量写明。
结果为:
-lambertw(-C1*exp(-x-1))-x-1
其中:
Y=lambertw(X)表示函数关系Y*exp(Y)=X。
例2:
求解常微分方程
的MATLAB程序为:
Y2=dsolve('y*D2y-Dy^2=0’,’x’)
结果为:
Y2=[exp((x+C2)/C1)]
[C2]
我们看到有两个解,其中一个是常数。
例3:
求常微分方程组
通解的MATLAB程序为:
[X,Y]=dsolve('Dx+5*x+y=exp(t),Dy-x-3*y=exp(2*t)','t')
例4:
求常微分方程组
通解的MATLAB程序为:
[X,Y]=dsolve('Dx+2*x-Dy=10*cos(t),Dx+Dy+2*y=4*exp(-2*t)','x(0)=2','y(0)=0')
以上这些都是常微分方程的精确解法,也称为常微分方程的符号解。
但是,我们知道,有大量的常微分方程虽然从理论上讲,其解是存在的,但我们却无法求出其解析解,此时,我们需要寻求方程的数值解,在求常微分方程数值解方面,MATLAB具有丰富的函数,我们将其统称为solver,其一般格式为:
[T,Y]=solver(odefun,tspan,y0)
该函数表示在区间tspan=[t0,tf]上,用初始条件y0求解显式常微分方程
。
solver为命令ode45,ode23,ode113,ode15s,ode23s,ode23t,ode23tb之一,这些命令各有特点。
我们列表说明如下:
求解器
ODE类型
特点
说明
ode45
非刚性
一步算法,4,5阶Runge-Kutta
方法累积截断误差
大部分场合的首选算法
ode23
非刚性
一步算法,2,3阶Runge-Kutta
方法累积截断误差
使用于精度较低的情形
ode113
非刚性
多步法,Adams算法,高低精度均可达到
计算时间比ode45短
ode23t
适度刚性
采用梯形算法
适度刚性情形
ode15s
刚性
多步法,Gear’s反向
数值积分,精度中等
若ode45失效时,
可尝试使用
ode23s
刚性
一步法,2阶Rosebrock算法,
低精度。
当精度较低时,
计算时间比ode15s短
odefun为显式常微分方程
中的
tspan为求解区间,要获得问题在其他指定点
上的解,则令
(要求
单调),
y0初始条件。
例5:
求解常微分方程
,
,
的MATLAB程序如下:
fun=inline('-2*y+2*x*x+2*x');[x,y]=ode23(fun,[0,0.5],1)
结果为:
x=0,0.0400,0.0900,0.1400,0.1900,0.2400,0.2900,0.3400,0.3900,0.4400,0.4900,0.5000
y=1.0000,0.9247,0.8434,0.7754,0.7199,0.6764,0.6440,0.6222,0.6105,0.6084,0.6154,0.6179
例6:
求解常微分方程
的解,并画出解的图形。
分析:
这是一个二阶非线性方程,用现成的方法均不能求解,但我们可以通过下面的变换,将二阶方程化为一阶方程组,即可求解。
令:
,
,
,则得到:
接着,编写vdp.m如下:
functionfy=vdp(t,x)
fy=[x
(2);7*(1-x
(1)^2)*x
(2)-x
(1)];
再编写m文件sy12_6.m如下:
y0=[1;0]
[t,x]=ode45(@vdp,[0,40],y0);
y=x(:
1);dy=x(:
2);
plot(t,y,t,dy)
已知一阶常微分方程,
g+0.047*du/dt+u/6.7=(13-u)/18
其中,当sin(10/pi*t)>=0时,g=1.18sin(10/pi*t)
当sin(10/pi*t)<0时,g=0
u的初值为u(0)=0,求t>=0时的解。
画出图像,并求u的极大值。
曾经将g写成,0.59sin(10/pi*t)+abs(0.59*sin(10/pi*t)),用dsolve求解,可惜出错了。
请大家帮忙,谢谢!
(用其他函数解也可以。
)
程序如下
fun=inline(['((13-u)/18-(sin(10*t/pi)>0)*',...
'1.18*sin(10*t/pi)-u/6.7)/0.047'],'t','u');
[t,u]=ode45(fun,[0,10],[0]);
plot(t,u)
说明g这样表示的:
gt=(sin(10*t/pi)>0)*1.18*sin(10*t/pi);
matlab一阶微分方程的解法:
Examples:
dsolve('Dx=-a*x')returns
ans=exp(-a*t)*C1
x=dsolve('Dx=-a*x','x(0)=1','s')returns
x=exp(-a*s)
y=dsolve('(Dy)^2+y^2=1','y(0)=0')returns
y=
[sin(t)]
[-sin(t)]
S=dsolve('Df=f+g','Dg=-f+g','f(0)=1','g(0)=2')
returnsastructureSwithfields
S.f=exp(t)*cos(t)+2*exp(t)*sin(t)
S.g=-exp(t)*sin(t)+2*exp(t)*cos(t)
Y=dsolve('Dy=y^2*(1-y)')
Warning:
Explicitsolutioncouldnotbefound;implicitsolutionreturned.
Y=
t+1/y-log(y)+log(-1+y)+C1=0
dsolve('Df=f+sin(t)','f(pi/2)=0')
dsolve('D2y=-a^2*y','y(0)=1,Dy(pi/a)=0')
S=dsolve('Dx=y','Dy=-x','x(0)=0','y(0)=1')
S=dsolve('Du=v,Dv=w,Dw=-u','u(0)=0,v(0)=0,w(0)=1')
w=dsolve('D3w=-w','w(0)=1,Dw(0)=0,D2w(0)=0')
y=dsolve('D2y=sin(y)');pretty(y)
12/22/2006
12.6一階微分方程之解法
12.6一階微分方程之解法
連續性函數之微分可以利用多種指令為之,而基本之細段微分法則是最早使用的概念,亦可利用MATLAB運算指令直接求解。
尤拉法(Eulermethod)
尤拉法為數值積分法最基本之一種。
它是利用時間細段微分,利用初值及一階微分方程進行求解。
以下列方程式為例:
y'=dy/dt=f(t,y)
式中t為時間,y為因變數。
y',dy/dt均為其一階微分的表示式。
m為t之函數。
根據微分之定義,可以改寫如下:
y'=dy/dt=limΔ-0[y(t+Δt)-y(t)]/dt=f(t,y)
展開,可以求得經過Δt之時間後,y(t+Δt)之結果為:
y(t+Δt)=y(t)+Δtf(t,y)
若設Δt為一均勻之區段時間,或稱為步進值,則可以給予y(t)一個初值,然後利用增加步進值之個數逐步累積至最終的區間,最後得到y之答案值。
其求解之型式如下:
y(tk+1)=y(tk)+Δtf(tk,y(tk))
故只要給定時間範圍及f(t0,y0)等初值,即可利用迴圈的方式得到最終答案。
其中,Δt設定值之大小則會影響答案之準確度,Δt值愈大,準確度會差,Δt值愈小則運算時間加長。
例:
設f(t,y)=-10y,初值t=0,y0=2。
且y'=-10y之真實解為y(t)=2e-10t。
設Δt=0.02,以此利用MATLAB求解。
functiony=euler(y0,time,delta)
%UsingEulermethodtosolvedifferentialeqs.
%Inputs:
%y0:
initialvalue
%time:
timelimit
%delta:
stepsize
%Example:
y=euler(2,0.5,0.02)
r=-10;k=0;
t=0:
delta:
time;y=zeros(size(t));
y
(1)=y0;
fori=2:
length(t),y(i)=y(i-1)+r*y(i-1)*delta;end
y_true=2*exp(-
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 例子