matlab解微分方程
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
应用差商 近似导数
dx f (t0 , x(t0 )) dt tt0 x(t0 h) x(t0 )
t0 h t0
x1 x0 hf (t0 , x0 ) R0
t0 ,
x(t0 )
舍入误差
•提高数值解精度的方法之一： 减小步长h的值；
•
•
function [outx,outy]=MyEuler(fun,x0,xt,y0,PointNum) % fun 表示f(x,y); x0,xt:自变量的初值和终值; %y0:函数在x0处的值,其可以为向量形式; %PointNum表示自变量在[x0,xt]上取的点数
if nargin==4 | PointNum<=0 %PointNum 默认值为100
PointNum=100; end if nargin==3 % y0默认值为0
故只有部分非线性微分方程可直接求出其解析解。
6.2Βιβλιοθήκη Baidu微分方程问题的数值解法
5.2.1 算法概述
dx
dt
f
(t, x),
x(a) x0.
a t b
• 微分方程求解的误差与步长问题：
a t0时刻系统状态向量的值为x(t0 ) 数值解法：即求微分方程初值问题的解在若干点
a t0 t1 t2 L tn b 的近似值y（k k 1, 2,..., n）的方法。 hk tk1 tk : 步长（一般取等步长h）
>> syms t x; x=dsolve('Dx=x*(1-x^2)+1')
Warning: Explicit solution could not be found; implicit solution returned. > In D:\MATLAB6p5\toolbox\symbolic\dsolve.m at line 292 x= t-Int(1/(a-a^3+1),a=``..x)+C1=0
数学描述：
d n y(t)
d n1 y(t)
d n2 y(t)
dy (t )
dt n a1 dt n1 a2 dt n2 L an1 dt an y(t) f (t)
特征方程（多项式代数方程）：
sn a1sn1 a2sn2 L an1s an 0
－该代数方程的根称为原常系数方程的特征根 －根据特征根的情况可求得原方程的解析解
• 例：求解
>> [x,y]=dsolve('D2x+2*Dx=x+2*y-exp(-t)', … 'Dy=4*x+3*y+4*exp(-t)')
• 例：
>> syms t x >> x=dsolve('Dx=x*(1-x^2)') x= [ 1/(1+exp(-2*t)*C1)^(1/2)] [ -1/(1+exp(-2*t)*C1)^(1/2)]
如果用推导的方法求Ci的值，每个系数的解析解至少要写出10 数行,故可采用有理式近似 的方式表示.
>> vpa(y,10) %有理近似值 ans = 1.196361839*exp(-5.*t)+.4166666667-
.4785447354*sin(t)*cos(t)*exp(-5.*t)-.4519262218e1*cos(2.*t)*exp(-5.*t)-2.392723677*cos(t)^2*exp(5.*t)+.2259631109*sin(2.*t)*exp(-5.*t)-473690.0893*exp(3.*t)+31319.63786*exp(-2.*t)-219.1293619*exp(1.*t)+442590.9059*exp(-4.*t)
分别处理系数，如：
>> [n,d]=rat(double(vpa(-445/26*cos(1)-51/13*sin(1)-69/2)))
ans = -8704
185 % rat()最接近有理数的分数
判断误差： >> vpa(-445/26*cos(sym(1))-51/13*sin(1)-69/2+8704/185) ans = .114731975864790922564144636e-4
• 格式： y=dsolve(f1, f2, …, fm)
• 格式：指明自变量 y=dsolve(f1, f2, …, fm ,’x’)
fi 即可以描述微分方程，又可描述初始条件 或边界条件。如：
描述微分方程时 y ( 4 ) (t ) 7 D 4 y 7
描述条件时
&y&(2) 3 D2 y(2) 3
例：
>> syms t; u=exp(-5*t)*cos(2*t+1)+5; >> uu=5*diff(u,t,2)+4*diff(u,t)+2*u uu = 87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10
>> syms t y; >> y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',...
'87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10'])
>> y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',...
'87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1) +10'], 'y(0)=3', 'Dy(0)=2', 'D2y(0)=0', 'D3y(0)=0')
第6章 微分方程问题的解法
• 常系数线性微分方程的解析解方法 • 常微分方程问题的数值解法
–微分方程问题算法概述 –四阶定步长 Runge-Kutta算法及 MATLAB 实现 –一阶微分方程组的数值解 –微分方程转换
• 特殊微分方程的数值解 • 边值问题的计算机求解 • 偏微分方程的解
6.1 常系数线性微分方程的解析解方法
>> y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',... '87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1) + ... 10'],'y(0)=1/2','Dy(pi)=1','D2y(2*pi)=0','Dy(2*pi)=1/5')
dx f (t0 , x(t0 )) dt tt0 x(t0 h) x(t0 )
t0 h t0
x1 x0 hf (t0 , x0 ) R0
t0 ,
x(t0 )
舍入误差
•提高数值解精度的方法之一： 减小步长h的值；
•
•
function [outx,outy]=MyEuler(fun,x0,xt,y0,PointNum) % fun 表示f(x,y); x0,xt:自变量的初值和终值; %y0:函数在x0处的值,其可以为向量形式; %PointNum表示自变量在[x0,xt]上取的点数
if nargin==4 | PointNum<=0 %PointNum 默认值为100
PointNum=100; end if nargin==3 % y0默认值为0
故只有部分非线性微分方程可直接求出其解析解。
6.2Βιβλιοθήκη Baidu微分方程问题的数值解法
5.2.1 算法概述
dx
dt
f
(t, x),
x(a) x0.
a t b
• 微分方程求解的误差与步长问题：
a t0时刻系统状态向量的值为x(t0 ) 数值解法：即求微分方程初值问题的解在若干点
a t0 t1 t2 L tn b 的近似值y（k k 1, 2,..., n）的方法。 hk tk1 tk : 步长（一般取等步长h）
>> syms t x; x=dsolve('Dx=x*(1-x^2)+1')
Warning: Explicit solution could not be found; implicit solution returned. > In D:\MATLAB6p5\toolbox\symbolic\dsolve.m at line 292 x= t-Int(1/(a-a^3+1),a=``..x)+C1=0
数学描述：
d n y(t)
d n1 y(t)
d n2 y(t)
dy (t )
dt n a1 dt n1 a2 dt n2 L an1 dt an y(t) f (t)
特征方程（多项式代数方程）：
sn a1sn1 a2sn2 L an1s an 0
－该代数方程的根称为原常系数方程的特征根 －根据特征根的情况可求得原方程的解析解
• 例：求解
>> [x,y]=dsolve('D2x+2*Dx=x+2*y-exp(-t)', … 'Dy=4*x+3*y+4*exp(-t)')
• 例：
>> syms t x >> x=dsolve('Dx=x*(1-x^2)') x= [ 1/(1+exp(-2*t)*C1)^(1/2)] [ -1/(1+exp(-2*t)*C1)^(1/2)]
如果用推导的方法求Ci的值，每个系数的解析解至少要写出10 数行,故可采用有理式近似 的方式表示.
>> vpa(y,10) %有理近似值 ans = 1.196361839*exp(-5.*t)+.4166666667-
.4785447354*sin(t)*cos(t)*exp(-5.*t)-.4519262218e1*cos(2.*t)*exp(-5.*t)-2.392723677*cos(t)^2*exp(5.*t)+.2259631109*sin(2.*t)*exp(-5.*t)-473690.0893*exp(3.*t)+31319.63786*exp(-2.*t)-219.1293619*exp(1.*t)+442590.9059*exp(-4.*t)
分别处理系数，如：
>> [n,d]=rat(double(vpa(-445/26*cos(1)-51/13*sin(1)-69/2)))
ans = -8704
185 % rat()最接近有理数的分数
判断误差： >> vpa(-445/26*cos(sym(1))-51/13*sin(1)-69/2+8704/185) ans = .114731975864790922564144636e-4
• 格式： y=dsolve(f1, f2, …, fm)
• 格式：指明自变量 y=dsolve(f1, f2, …, fm ,’x’)
fi 即可以描述微分方程，又可描述初始条件 或边界条件。如：
描述微分方程时 y ( 4 ) (t ) 7 D 4 y 7
描述条件时
&y&(2) 3 D2 y(2) 3
例：
>> syms t; u=exp(-5*t)*cos(2*t+1)+5; >> uu=5*diff(u,t,2)+4*diff(u,t)+2*u uu = 87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10
>> syms t y; >> y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',...
'87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10'])
>> y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',...
'87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1) +10'], 'y(0)=3', 'Dy(0)=2', 'D2y(0)=0', 'D3y(0)=0')
第6章 微分方程问题的解法
• 常系数线性微分方程的解析解方法 • 常微分方程问题的数值解法
–微分方程问题算法概述 –四阶定步长 Runge-Kutta算法及 MATLAB 实现 –一阶微分方程组的数值解 –微分方程转换
• 特殊微分方程的数值解 • 边值问题的计算机求解 • 偏微分方程的解
6.1 常系数线性微分方程的解析解方法
>> y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',... '87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1) + ... 10'],'y(0)=1/2','Dy(pi)=1','D2y(2*pi)=0','Dy(2*pi)=1/5')