maple 解代数方程组得多项式_Maple笔记2--常微分方程求解
来源:网络论坛转载(VB资料库)
常微分方程求解
微分方程求解是数学研究与应用的一个重点和难点. Maple能够显式或隐式地解析地求解许多微分方程求解.
在常微分方程求解器dsolve中使用了一些传统的技术例如laplace变换和积分因子法等,
函数pdesolve则使用诸如特征根法等经典方法求解偏微分方程. 此外, Maple还提供了可作摄动解的所有工具,
例如Poincare-Lindstedt法和高阶多重尺度法.
帮助处理常微分方程(组)的各类函数存于Detools软件包中, 函数种类主要有:可视化类的函数, 处理宠加莱动态系统的函数,
调整微分方程的函数, 处理积分因子、李对称法和常微分方程分类的函数, 微分算子的函数, 利用可积性与微分消去的方法简化微分方程的函数,
以及构造封闭解的函数等.
更重要的是其提供的强大的图形绘制命令Deplot能够帮助我们解决一些较为复杂的问题.
2.1 常微分方程的解析解
求解常微分方程最简单的方法是利用求解函数dsolve. 命令格式为:
dsolve(ODE);
dsolve(ODE, y(x),
extra_args);
dsolve({ODE, ICs}, y(x),
extra_args);
dsolve({sysODE, ICs}, {funcs}, extra_args);
其中, ODE—常微分方程, y(x)—单变量的任意变量函数, Ics—初始条件, {sysODE}—ODE方程组的集合,
{funcs}—变量函数的集合,
extra_args—依赖于要求解的问题类型.
例如, 对于一阶常微分方程 可用dsolve直接求得解析解:
> ODE:=x*diff(y(x),x)=y(x)*ln(x*y(x))-y(x);
> dsolve(ODE,y(x));
可以看出, dsolve的第一个参数是待求的微分方程, 第二个参数是未知函数. 需要注意的是, 无论在方程中还是作为第二个参数,
未知函数必须用函数的形式给出(即:必须加括号, 并在其中明确自变量), 这一规定是必须的,
否则Maple将无法区分方程中的函数、自变量和参变量, 这一点和我们平时的书写习惯不一致. 为了使其与我们的习惯一致,
可用alias将函数用别称表示:
> alias(y=y(x));
> ODE:=x*diff(y,x)=y*ln(x*y)-y;
> dsolve(ODE,y);
函数dsolve给出的是微分方程的通解,
其中的任意常数是用下划线起始的内部变量表示的.
在Maple中, 微分方程的解是很容易验证的,
只需要将解代入到原方程并化简就可以了.
> subs(%,ODE);
> assume(x,real): assume(_C1,real):
> simplify(%);
> evalb(%);
evalb函数的目的是对一个包含关系型运算符的表达式使用三值逻辑系统求值, 返回的值是true, false和FAIL.
如果无法求值, 则返回一个未求值的表达式. 通常包含关系型运算符“=, <>, ,
>=”的表达式在Maple中看作是代数方程或者不等式. 然而,
作为参数传递给evalb或者出现在if或while语句的逻辑表达式中时, 它们会被求值为true或false. 值得注意的是,
evalb不化简表达式, 因此在使用evalb之前应将表达式化简, 否则可能会出错.
再看下面常微分方程的求解:
> alias(y=y(x)):
> ODE:=diff(y,x)=sqrt(y^2+1);
> dsolve(ODE,y);
函数dsolve对于求解含有未知参变量的常微分方程也完全可以胜任:
> alias(y=y(x)):
> ODE:=diff(y,x)=-y/sqrt(a^2-y^2);
> sol:=dsolve(ODE,y);
由此可见, 对于不能表示成显式结果的微分方程解, Maple尽可能将结果表示成隐式解. 另外, 对于平凡解y=0常常忽略,
这一点应该引起注意.
dsolve对于求解微分方程初值问题也十分方便的:
> ODE:=diff(u(t),t$2)+omega^2*u(t)=0;
> dsolve({ODE,u(0)=u0,D(u)(0)=v0},u(t));
2.2 利用积分变换求解微分方程
对于特殊的微分方程, 我们还可以指定dsolve利用积分变换方法求解,
只需要在dsolve中加入可选参数method=transform即可. 其中transform是积分变换,
可以是laplace、fourier、fouriercos或者fouriersin变换.
作为例子, 我们来看一个具有阻尼的振子在阶跃冲击(Heaviside函数)下的响应:
>
ODE:=diff(u(t),t$2)+2*d*omega*diff(u(t),t)+omega^2*u(t)=Heaviside(t);
> initvals:=(u(0)=u[0],D(u)(0)=v[0]);
>
solution:=dsolve({ODE,initvals},u(t),method=laplace);
Maple给出了问题的通解, 但没有区分自由振动(d=0)、欠阻尼(01)的情况.
下面加以区分求解:
> assume(omega>0):
> simplify(subs(d=0,solution));
> K:=subs(d=1/5,u[0]=1,v[0]=1,solution);
> with(plots):
>
plot3d(rhs(%%),omega=2/3..4/3,t=0..20,style=hidden,orientation=[-30,45],axes=framed);
对于d=1的情况, 可可用下式获得结果:
> limit(rhs(solution),d=1);
再如下例:
>
diff(u(t),t$2)+3*diff(u(t),t)+2*u(t)=exp(-abs(t));
> dsolve(%,u(t),method=fourier);
2.3 常微分方程组的求解
函数dsolve不仅可以用来求解单个常微分方程, 也可以求解联立的常微分方程组. 特别是对于线性微分方程组,
由于数学上具有成熟的理论, Maple的求解也是得心应手. 其命令格式为:
dsolve( {eqn1, eqn2, …, ini_conds}, {vars});
其中, ini_conds是初始条件.
>
eqn1:={diff(x(t),t)=x(t)+y(t),diff(y(t),t)=y(t)-x(t)};
> dsolve(eqn1,{x(t),y(t)});
> eqn2:=2*diff(x(t),t$2)+2*x(t)+y(t)=2*t;
> eqn3:=diff(y(t),t$2)+2*x(t)+y(t)=t^2+1;
> dsolve({eqn2, qn3, x(0)=0, D(x)(0)=1, y(0)=0, D(y)(0)=0},
{x(t),y(t)} );
2.4 常微分方程的级数解法
1) 泰勒级数解法
当一个常微分方程的解析解难以求得时, 可以用Maple求得方程解的级数近似, 这在大多数情况下是一种非常好的方法.
级数解法是一种半解析半数值的方法.
泰勒级数法的使用命令为:
dsolve({ODE,Ics}, y(x), 'series'); 或dsolve({ODE,Ics}, y(x),
'type=series');
下面求解物理摆的大幅振动方程: , 其中l是摆长, 是摆角,
g是重力加速度.
> ODE:=l*diff(theta(t),t$2)=-g*sin(theta(t));
> initvals:=theta(0)=0,D(theta)(0)=v[0]/l;
>
sol:=dsolve({ODE,initvals},theta(t),type=series);
> Order:=11:
>
sol:=dsolve({ODE,initvals},theta(t),type=series);
2) 幂级数解法
对于一个符号代数系统来说, 幂级数是必不可少的微分方程求解工具. 幂级数求解函数powsolve存于工具包powseries中.
但是, 这一求解函数的使用范围很有限,
它只可以用来求解多项式系数的线性常微分方程或方程组,其求解命令为:powseries[function]
(prep)或直接载入软件包后用function(prep),
prep为求解的线性微分方程及其初值.
例:求解:
>
ODE:=x*diff(y(x),x$2)+diff(y(x),x)+4*x^2*y(x)=0;
> dsolve(ODE,y(x));
> initvals:=y(0)=y0,D(y)(0)=0;
> with(powseries):
> sol:=powsolve({ODE,initvals});
> tpsform(sol,x,16);
也可以用powsolve给出的函数直接获得用递归形式定义的幂级数系数, 不过参数必须用_k,
这是powsolve使用的临时变量.
> sol(_k);
例:求解一维谐振子的解:
> alias(y=y(x)):
> ODE:=diff(y,x$2)+(epsilon-x^2)*y=0;
> H:=powsolve(ODE);
> tpsform(H,x,8);
> H(_k);
2.5 常微分方程的数值解法
在对微分方程的解析解失效后, 可以求助于数值方法求解微分方程. 数值求解的好处是只要微分方程的条件足够多时一般都可求得结果,
然而所得结果是否正确则必须依赖相关数学基础加以判断.
调用函数dsolve求常微分方程初值问题的数值解时需加入参数type=numeric.
另一方面, 常微分方程初值问题数值求解还可以选择算法, 加入参数“method=方法参数”即可,
方法参数主要有:
rkf45:4~5阶变步长Runge-Kutta-Fehlberg法
dverk78:7~8阶变步长Runge-Kutta-Fehlberg法
classical:经典方法, 包括向前欧拉法, 改进欧拉法, 2、3、4阶龙格库塔法,
Sdams-Bashford方法等
gear:吉尔单步法
mgear:吉尔多步法
2.5.1变步长龙格库塔法
下面用4~5阶Runge-Kutta-Fehlberg法求解van der Pol方程:
>
ODE:=diff(y(t),t$2)-(1-y(t)^2)*diff(y(t),t)+y(t)=0;
> initvals:=y(0)=0,D(y)(0)=-0.1;
> F:=dsolve({ODE,initvals},y(t),type=numeric);
此时, 函数返回的是一个函数,
可以在给定的数值点上对它求值:
> F(0);
> F(1);
可以看到,
F给出的是一个包括t、y(t)、D(y)(t)在内的有序表, 它对于每一个时间点可以给出一组数值表达式. 有序表的每一项是一个等式,
可对其作图描述.
> plot('rhs(F(t)[2])', t=0..15, title="solution of the Van de
Pol's Equation");
> plots[odeplot](F,[t,y(t)],0..15,title="solution of the Van de
Pol's Equation");
2.5.2吉尔法求解刚性方程
在科学和工程计算中, 常常会遇到这样一类常微分方程问题, 它可以表示成方程组: , 称其为刚性方程, 其解的分量数量相差很大,
分量的变化速度也相差很大. 如果用常规方法求解, 为了使变量有足够高的精度, 必须取很小的步长, 而为了使慢变分量达到近似的稳态解,
则需要很长的时间, 这样用小步长大时间跨度的计算, 必定造成庞大的计算量, 而且会使误差不断积累.
吉尔法是专门用来求解刚性方程的一种数值方法.
>
ODE:=diff(u(t),t)=-2000*u(t)+999.75*v(t)+1000.25,diff(v(t),t)=u(t)-v(t);
> initvals:=u(0)=0,v(0)=-2;
>
ansl:=dsolve({ODE,initvals},{u(t),v(t)},type=numeric,method=gear);
> ansl(10,0);
> p1:=plots[odeplot]
(ansl,[t,u(t)],0..20,color=red):
p2:=plots[odeplot]
(ansl,[t,v(t)],0..20,color=blue):
plots[display] ({p1,p2}, title="Solution of a stiff
equation");
2.5.3 经典数值方法
Maple中常微分方程数值解法中有一类被称作是“经典”(classical)方法. 当然, 称其为经典方法不是因为它们常用或是精度高,
而是因为它们的形式简单, 经常被用于计算方法课上的教学内容. 它们是一些常见的固定步长方法,
在dsolve中用参数method=classical[方法名称], 如果不特别指出, 将默认采用向前欧拉法.
主要有:
foreuler:向前欧拉法(默认)
hunform:Heun公式法(梯形方法, 改进欧拉法)
imply:改进多项式法
rk2:二阶龙格库塔法
rk3:三阶龙格库塔法
rk4:四阶龙格库塔法
adambash:Adams-Bashford方法(预测法)
abmoulton:Adams-Bashford-Moulton方法(预测法)
下面给出微分方程数值方法的参数表:
参数名
参数类型
参数用途
参数用法
initial
浮点数的一维数组
指定初值向量
number
正整数
指定向量个数
output
'procedurelist'(默认)
或'listprocedure'
指定生成单个函数或多个函数的有序表
Procedurelis:单个函数,返回有序表
Listprocedure:函数的有序表
procedure
子程序名
用子程序形式指定第一尖常微分方程组的右边部分
参数1:未知函数的个数
参数2:自变量
参数3:函数向量
参数4:导函数向量
start
浮点数
自变量起始值
startinit
布尔量(默认FALSE)
指定数值积分是否总是从起始值开始
对dverk78不适用
value
浮点数向量(一维数组)
指定需要输出函数值的自变量数值点
如果给定,结果是一个的矩阵.元素[1,1]是一个向量,含自变量名和函数名称;元素[2,1]是一个数值矩阵,其中第一列value的输入相同,其他列中是相应的函数值
另外, 还有一些特殊的附加参数:
maxfun:整数类型, 用于最大的函数值数量, 默认值50000, 为负数时表示无限制
corrections:正整数类型, 指定每步修正值数量, 在abmoulton中使用, 建议值≤4
stepsize:浮点数值, 指定步长
下面看一个简单的例子:
> ODE:=diff(y(x),x)=y(x)-2*x/y(x);
> initvals:=y(0)=1;
>
sol1:=dsolve({ODE,initvals},y(x),numeric,method=classical,stepsize=0.1,start=0);
而其解析解为:
> sol2:=dsolve({diff(y(x),x)=y(x)-2*x/y(x), y(0)=1},
y(x));
将两者图形同时绘制在同一坐标系中比较, 可以发现,
在经过一段时间后, 欧拉法的数值结果会产生较大误差.
> plot({rhs(sol2),'rhs(sol1(x)[2])'},x=0..2);
求解微方程, 无论使用什么方法或者加入什么选项,
求解完成后必须利用相关数学知识进行逻辑判断, 绝对不对简单迷信Maple给出的结果, 否则很有可能得到一个对于方程本身也许还看得过去,
但在数学或者物理意义上不合理的解.
2.6摄动法求解常微分方程
由于微分方程求解的复杂性, 一般微分方程常常不能求得精确解析解, 需要借助其它方法求得近似解或数值解, 或者两种方法兼而有之.
摄动法是重要的近似求解方法.
摄动法又称小参数法, 它处理含小参数 的系统, 一般当 =0时可求得解x0. 于是可把原系统的解展成 的幂级数 ,
若这个级数当 0时一致收敛,则称正则摄动, 否则称奇异摄动.
摄动法的种类繁多, 最有代表性的是庞加莱—林斯泰特(Poicare-Lindstedt)法, 在此, 我们以该方法求解van der
Pol方程:
当 =0时该方程退化为数学单摆的常微分方程, 当 =1时为3.5讨论的情况, 对任意 , 该微分方程拥有一个渐进稳定的周期解,
称为极限环.
由于van der Pol方程中没有显式的时间依赖项, 不失一般性, 设初值为y(0)=0. 在庞加莱—林斯泰特法中,
时间通过变换拉伸:
, 其中
对于 , van der Pol方程变为:
restart:
diff(y(t),t$2)-epsilon*(1-y(t)^2)*diff(y(t),t)+y(t)=0;
>
ODE:=DEtools[Dchangevar]({t=tau/omega,y(t)=y(tau)},%,t,tau);
> e_order:=6:
> macro(e=epsilon,t=tau):
> alias(seq(y=eta(tau),i=0..e_order)):
> e:=()->e:
> for i from 0 to e_order do
eta:=t->eta(t)
od:
> omega:=1+sum('w*e^i','i'=1..e_order);
> y:=sum('eta*e^i','i'=0..e_order);
> deqn:=simplify(collect(ODE,e),{e^(e_order+1)=0}):
> for i from 0 to e_order do
ode:=coeff(lhs(deqn),e,i)=0
od:
> ode[0];
> ode[1];
> ode[2];
>
dsolve({ode[0],eta[0](0)=0,D(eta[0])(0)=C[1]},eta[0](t));
> eta[0]:=unapply(rhs(%),t);
> ode[1];
> map(combine,ode[1],'trig');
> ode[1]:=map(collect,%,[sin(t),cos(t)]);
>
dsolve({ode[1],eta[1](0)=0,D(eta[1])(0)=C[2]},eta[1](t),method=laplace);
> map(collect,%,[sin(t),cos(t),t]);
>
solve({coeff(lhs(ode[1]),sin(t))=0,coeff(lhs(ode[1]),cos(t))=0});
> w[1]:=0:C[1]:=-2:
> ode[1];
>
dsolve({ode[1],eta[1](0)=0,D(eta[1])(0)=C[2]},eta[1](t),method=laplace);
> eta[1]:=unapply(rhs(%),tau);
> map(combine,ode[2],'trig'):
>
ode[2]:=map(collect,%,[sin(t),sin(3*t),cos(t),cos(3*t)]);
>
solve({coeff(lhs(ode[2]),sin(t))=0,coeff(lhs(ode[2]),cos(t))=0});
> assign(%):
>
dsolve({ode[2],eta[2](0)=0,D(eta[2])(0)=C[3]},eta[2](t),method=laplace):
> eta[2]:=unapply(rhs(%),t);
> for i from 0 to e_order do
map(combine,ode,'trig'):
ode:=map(collect,%,[seq(sin((2*j+1)*t),j=0..i),seq(cos((2*j+1)*t),j=0..i)]):
solve({coeff(lhs(ode),sin(t))=0,coeff(lhs(ode),cos(t))=0}):
assign(%):
dsolve({ode,eta(0)=0,D(eta)(0)=C[i+1]},eta(t),method=laplace):
collect(%,[seq(sin((2*j+1)*t),j=0..i),seq(cos((2*j+1)*t),j=0..i)]):
eta:=unapply(rhs(%),t);
od:
> omega;
> y(t):
> y:=unapply(simplify(y(t),{e^e_order=0}),t):
> e:=1:y(t);
> plot(y(t),t=0..15);
在该问题的求解过程中, 前半部分我们按照交互式命令方式输入, 也就是把数学逻辑推理的过程“翻译”成Maple函数, 而在后半部分,
则采用程序设计方式书写了数学推导过程, 这是应用Maple解决实际问题的两种方式. 前一种方法只需了解Maple函数即可应用,
而后一种程序设计方式则需掌握Maple程序设计语言. 但是, 不论是那一种方式, 数学基础总是最重要的.
maple 解代数方程组得多项式_Maple笔记2--常微分方程求解相关推荐
- maple 解代数方程组得多项式_Maple与数学实验
Maple与数学实验 出版时间:2013年版 内容简介 数学实验是利用数学软件借助于计算机来处理数学问题的一门学科.<Maple与数学实验>介绍了Maple软件在符号运算.数值计算以及绘图 ...
- maple 解代数方程组得多项式_利用修正影射法求组合KdV方程新的精确解
1引言非线性科学研究的一个重要方面就是讨论孤立子的性质.相互作用及其随时间运动演化的特点,因此非线性演化方程的求解越来越显得具有理论和实际意义.组合KdV方程是KdV和mKdV方程的复合,既包含有非线 ...
- matlab中欠定方程组超定方程组_MATLAB解代数方程组一些函数用法1
1.solve函数用法 solve('函数方程组')---解方程 ezplot('函数方程组',[x1 x2 y1 y2])---画函数的方程 root(f,x,k)--f表达式的k阶开根,x是变量. ...
- MATLAB 解代数方程组
今天用matlab给了我三重惊喜,简直打开了新世界的大门: 1.虽然知道matlab有内置的符号工具箱,但以往用的很少,直到今天,需要求解一个方程组,方程本身到不是多么复杂,只不过变量众多,非常的恶心 ...
- matlab求解代数方程组,matlab求解代数方程组.doc
matlab求解代数方程组.doc 1第三讲Matlab求解代数方程组理论介绍:直接法+迭代法,简单介绍相关知识和应用条件及注意事项软件求解:各种求解程序讨论如下表示含有个未知数.由个方程构成的线性方 ...
- 高等数值计算方法学习笔记第6章【解线性代数方程组的迭代方法(高维稀疏矩阵)】
高等数值计算方法学习笔记第6章[解线性代数方程组的迭代方法(高维稀疏矩阵)] 一.引言 1.例题(说明迭代法的收敛性研究的重要性) 2.定义(迭代法,迭代法收敛)&解误差 二.基本迭代法 1. ...
- 如何用matlab解异或方程,Matlab-6:解非线性方程组newton迭代法
函数文件: function x=newton_Iterative_method(f,n,Initial) x0=Initial; tol=1e-11; x1=x0-Jacobian(f,n,x0)\ ...
- Maple笔记2--常微分方程求解
转需看 原文地址:Maple笔记2--常微分方程求解作者:Lionel 来源:网络论坛转载(VB资料库) 常微分方程求解 微分方程求解是数学研究与应用的一个重点和难点. Maple能够显式或隐式地解析 ...
- 第二十八讲 解非齐次线性方程组
一,关于二阶方程组x⃗′=Ax⃗{\vec{x}}'=A\vec{x}x′=Ax的理论(对n阶方程也成立): (假设A是常数矩阵) 定理A:x⃗′=Ax⃗{\vec{x}}'=A\vec{x}x′=A ...
最新文章
- 关于php无字母代码的研究
- JavaScript 技术篇-简单的两行js代码获取password不可见密码实例演示,js获取密码输入框里的值
- 基于内容的图像检索 Database for Content-Based Image Retrieval
- 【NLP】从头开始学词向量的预训练
- 【Qt】2D绘图之窗口-视口转换
- phpmyadmin教程:使用phpmyadmin创建用户
- M. Monster Hunter(树形dp)
- 【APICloud系列|14】xcode下载地址
- 物资管理系统c语言课程设计,C语言实现仓库物资管理系统
- leetcode 51. N 皇后 思考分析
- 四张图,读懂 BIO、NIO、AIO、多路复用 IO 的区别
- DevOps实践教程 华为云 系列教程2021 合集
- 【Shell脚本学习7】Shell脚本学习指南分享
- 史上最全GIS相关软件(CAD、FME、Arcgis、ArcgisPro)
- 电脑显示器黑屏|显示器突然黑屏|显示器闪黑屏
- Linux(Ubuntu)菜单栏(工具栏)隐藏了,怎么显示出来
- LTP性能测试工具的使用详解
- 在线演绎3D图表如何操作
- 两分钟研究透idea中Git文件的颜色,绿红蓝白灰
- CMake入门教程【手册篇】CMake生成与编译项目
热门文章
- 《Discriminative Unsupervised Feature Learning with Exemplar Convolutional Neural Networks》阅读笔记
- pet-shop Dapp开发(下)
- 在树莓派上面 玩 仙剑奇侠传
- linux 拼图游戏,立体艺术拼图游戏
- 清风算法对seo不是打击而是好事
- 物联网毕业设计 单片机智能扫地机器人设计与实现
- 树莓派+Ardunio的魔方机器人
- 阿里云邮箱登录日志中有异地IP登录是怎么回事?该怎么办?
- 滴答顺风车怎么抢90%以上的订单_想来赚顺风车钱的补课内容都给你准备好了
- JavaScript 数组塌陷