MATLAB编程基础第6讲--插值、拟合与初值常微分方程的求解

matlab 数学建模

3.9常微分方程初值问题的数值解法 龙格-库塔法简介 龙格-库塔法的实现 基于龙格-库塔法，MATLAB提供了求常微分方程数值解 的函数，一般调用格式为： [t,y]=ode23('fname',tspan,y0) [t,y]=ode45('fname',tspan,y0) 其中fname是定义f(t,y)的函数文件名，该函数文件必须返回 一个列向量。tspan形式为[t0,tf],表示求解区间。y0是初始 状态列向量。t和y分别给出时间向量和相应的状态向量。8/12/2013 1

matlab 数学建模

一阶常微分方程的初值问题求解

y f x, y , a x b y x0 y0 y1 y2 y3 例3-39,求微分方程组 y2 y1 y3 在区间[0 y 0.51y y 1 2 312]内的数值解，且满足初始条件

y1 0 0 y2 0 1 y 0 1 32

8/12/2013

matlab 数学建模

例3-39 % 首先编写方程函数，名为rigid.m function dy = rigid(t,y) dy = zeros(3,1); % 生成3×1的矩阵 dy(1) = y(2) * y(3); % 第一个元素是 dy(2) = -y(1) * y(3); % 第二个元素是 dy(3) = -0.51 * y(1) * y(2); % 第三个元素是

8/12/2013

matlab 数学建模

例3-39(续前) tspan = [0 12]; y0 = [0 1 1]; % 在同一目录下，计算方程数值解。输入时间区间和初始 条件 [t,Y] = ode45('rigid',tspan,y0); % 采用ode45算法求解方程，options为默认值 plot(t,Y(:,1),'-',t,Y(:,2),'-.',t,Y(:,3),'.') % 绘制计算结果并标注，如图3-14 legend('Y1','Y2','Y3')8/12/2013 4

matlab 数学建模

高阶微分方程的初值问题求解F y, y , y , , y n , x 0高阶微分方程：

y 0 y0 , y 0 y0 , y n 1 0 y0n 1

把高阶微分方程转换为一阶方程组，初始条件也做相应的替换。通常令： n 1 上式改写为： 1 2 n

y y, y y , , y y y1 f1 x, Y y f x, Y Y 2 2 yn f n x, Y

y1 0 y0 y 0 y Y0 2 0 n 1 yn 0 y0 5

8/12/2013

matlab 数学建模

例3-40,求方程 y 2 y 2 y 4 y1 2 0 初始条件

y 0 1, y 0 1, y 0 0

在时间区间[0 12]内的数值解

把高阶微分转化低阶微分形式，令：

y1 y, y2 y , y3 y

y1 y2 y2 y3 原方程变为： y y 1 2 2 y 2 y 4 1 3 2 3

y1 0 1 初始条件： y2 0 1 y 0 0 38/12/2013

matlab 数学建模

例3-40 训练任务：请编制名为odet3.m方程函数 tspan = [0 12]; % 在同一目录下，计算方程数值解。给出时间区间和

初始 值 y0 = [-1 1 0]; [t,Y] = ode45('odet3',tspan,y0); % 采用ode45算法 plot(t,Y(:,1),'-',t,Y(:,2),'-.',t,Y(:,3),'.') % 绘制计算结果并标注，如图3-15 legend('Y1','Y2','Y3')8/12/2013 7

matlab 数学建模

% 首先编写方程函数，名为odet3.m function dy = odet3(t,y) dy = zeros(3,1); dy(1) = y(2); dy(2) = y(3); dy(3) = y(1).^(1/2)-2*y(3).^2*(y(2)-4);

8/12/2013

matlab 数学建模

例设有初值问题，试求其数值解，并与精确解相比较 (精确解为y(t)=)。(1) 建立函数文件funt.m。 function yp=funt(t,y) yp=(y^2-t-2)/4/(t+1); (2) 求解微分方程。 t0=0;tf=10; y0=2; [t,y]=ode23('funt',[t0,tf],y0); %求数值解 y1=sqrt(t+1)+1; %求精确解 t' y' y1' 8/12/2013 y为数值解，y1为精确值，显然两者近似。

matlab 数学建模

· 2-1)x+x=0 例：x+(x

· ·

为方便令x1=x,x2=x分别对x1,x2求一

·

阶导数，整理后写成一阶微分方程组形式· x1=x2 · x2=x2(1-x12)-x1

1. 建立m文件8/12/2013

2. 解微分方程

matlab 数学建模

建立m文件 function xdot=wf(t,x) xdot=zeros(2,1) xdot(1)=x(2) xdot(2)=x(2)*(1-x(1)^2)-x(1) 给定区间、初始值;求解微分方程 t0=0; tf=20; x0=[0 0.25]'; [t,x]=ode23('wf', t0, tf, x0) figure(1), plot(t,x); figure(2),plot(x(:,1),x(:,2)) 8/12/2013

matlab 数学建模

3

2

1

0

-1

-2

-3

0

5

10

15

20

25

30

3

2

1

0

-1

-2

8/12/2013 -3

12-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5

-2.5

matlab 数学建模

训练13x 4 y 7 z 12w 4 5 x 7 y 4 z 2w 3 x 8 z 5w 9 6 x 5 y 2 z 10w 8a=[3 4 -7 -12;5 -7 4 2;1 0 8 -5;-6 5 -2 10]; b=[4 -3 9 -8]'; c=a\b

8/12/2013

matlab 数学建模

训练2算出 x3+2x2+x+1=0 的根 function f=fesin(x) f=exp(-0.5*x).*sin(x+pi/6); [S,n]=quad('fesin',0,3*pi) S= 0.9008 n= 77

训练3 3 0

e sin( x pi / 6) dx

x 2

a=quad('exp(-x/2).*sin(x+pi/6)',0,3*pi)

a =

0.9008

8/12/2013

matlab 数学建模

数据的保存与获取 把matlab工作空间中一些有用的数据 长久保存下来的方法是生成mat数据文 件。 save —— 将工作空间中所有的变量 存到matlab.mat文件中。默认文件名8/12/2013 15

matlab 数学建模

(1)save filename variables 将变量列表variables所列出的变量保存到磁盘文 件filename中 Variables所表示的变量列表中，不能用逗号，各 个不同的变量之间只能用空格来分隔。 未列出variables时，表示将当前工作空间中所 有变量都保持到磁盘文件中。(保存现场信息备 用) 缺省的磁盘文件扩展名为“.mat”,数据存储格 式为二进制。可以使用“-”定义不同的存储格 式。

8/12/2013