matlab求微分方程的解-实验分析报告四

matlab求微分方程的解-实验报告四

————————————————————————————————作者：————————————————————————————————日期：

《matlab与数学实验》实验报告

实验序号：实验四日期：2015年 5 月25 日

班级

132132002

姓名彭婉婷学号1321320057 实验名称求微分方程的解

问题背景描述

实际应用问题通过数学建模所归纳而得到的方程，绝大多数都是微分方程，另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）．这就要求我们既要研究微分方程（组）的解析解法（精确解），更要研究微分方程（组）的数值解法（近似解）．

实验目的本实验将主要研究微分方程(组)的数值解法（近似解），重点介绍Euler 折线法．

实验原理

与数学模

型

MATLAB7.11.0

（要点）

2.

求微分方程

x e y y y x sin 5'2''=+-的通解． 3. 求微分方程组

⎪⎪⎩⎪⎪⎨

⎧=-+=++00y x dt dy y x dt dx

在初始条件0|,1|00====t t y x 下的特解，并画出解函数()y f x =

的图形．

4. 分别用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解)，求解区间为

[0,2]t ∈．利用画图来比较两种求解器之间的差异．

5. 用 Euler 折线法求解微分方程初值问题

⎪⎩

⎪⎨

⎧=-=1)0(,

12'32y y x y y

的数值解(步长h 取0.1)，求解范围为区间[0,2]．

选做：

6. 用四阶 Runge-Kutta 法求解微分方程初值问题

⎩⎨

⎧=-=1)0(,cos 'y x e y y x

的数值解(步长h 取0.1)，求解范围为区间[0,3]．迭代法

记录（含基本步骤、主要程序清单及异常情况记录等）程序:

clear

syms x y

y=dsolve('(x^2-1)*Dy+2*x*y=sin(x)','x')

答案:y =-(C2 + cos(x))/(x^2 - 1)

2.求微分方程x

e

y

y

y x sin

5

'

2

''=

+

-的通解．

程序:

clear

syms x y

y=dsolve('D2y-2*Dy+5*y=exp(x)*sin(x)','x') simplyify(x/y)weijiao

答案:

y =(exp(x)*sin(x))/6 - (sin(3*x)*exp(x))/8 + (sin(5*x)*exp(x))/24 + C4*cos(2*x)*exp(x) + sin(2*x)*exp(x)*(cos(x)/4 - cos(3*x)/12 + 1/6) + C5*sin(2*x)*exp(x)

3. 求微分方程组

⎪⎪⎩⎪⎪⎨

⎧=-+=++00y x dt dy y x dt dx

在初始条件0|,1|00====t t y x 下的特解，并画出解函数()y f x =

的图形．

程序: clear syms x y t

[x,y]=dsolve('Dx+x+y=0','Dy+x-y=0','x(0)=1','y(0)=0','t') ezplot(x,y,[0, 1])(t 的取值，t 是与x,y 相关的，如果不给范围，则会默认为一个较大的区间) simplify(x) simplify(y) 答案:

x =exp(2^(1/2)*t)/2 + 1/(2*exp(2^(1/2)*t)) -

(2^(1/2)*exp(2^(1/2)*t))/4 + 2^(1/2)/(4*exp(2^(1/2)*t)) y =2^(1/2)/(4*exp(2^(1/2)*t)) - (2^(1/2)*exp(2^(1/2)*t))/4 图形:

4. 分别用ode23、ode45 求上述第3 题中的微分方程初值问题的数值解(近似解)，求解区间为[0,2]

t ．利用画图来比较两种求解器之间的差异．先编写函数文件verderpol.m：

clear

function xprime=verderpol(t,x)

xprime=[-x(1)-x(2); x(2)-x(1)];

再编写命令文件

clear

y0=[1;0];

[t,x] = ode45('verderpol',[0,2],y0);

x1=x(:,1);x2=x(:,2);

plot(x1,x2,'b-')

hold on

y0=[1;0];

[t,x]=ode23('verderpol',[0,2],y0); x1=x(:,1);x2=x(:,2);

plot(x(:,1),x(:,2),'r-');

图形:

两种求解器之间的差异:

ode45大部分场合的首选算法ode23使用于精度较低的情形但在此题中并没有体现出差异。

5. 用 Euler 折线法求解微分方程初值问题

⎪⎩

⎪⎨

⎧=-=1)0(,

12'32y y x y y

的数值解(步长h 取0.1)，求解范围为区间[0,2]．

程序：clear f=sym('y+2*x/y^2'); a=0; b=2; h=0.4; n=(b-a)/h+1; x=0; y=1; szj=[x,y]; for i=1:n-1

y=y+h*subs(f,{'x','y'},{x,y}); x=x+h; szj=[szj;x,y]; end szj

plot(szj(:,1),szj(:,2)) 答案：szj =