第一篇为 基础概念 ,第二篇为 R-K法的具体实现方法。

(一)常微分方程的MATLAB求解

概要:

  1. 常微分方程的MATLAB求解分为解析解、数值解
  2. 解析解(只有少数微分方程组有解析解):dsolve函数
  3. 数值解:solver函数,MATLAB内部自带了许多求解器(各种数值求解的方法,如 R-K法)
  4. 常用函数文件构造微分方程组
  5. 数值求解时:高阶微分方程化成常微分方程组

一、解析解

dsolve函数

  • y=dslove('eq1' , 'eq2' , ... , 'cond1' , 'cond2' , ... , 'v')
    • y——输出
    • eq1、eq2——微分方程
    • cond1、cond2——初值条件,省略时表示求通解
    • v——自变量,省略时默认为t
  • 微分方程中用 D 表示 自变量的导数
    • Dy—— y′
    • D2y—— y″
    • D3y—— y‴
  • 若找不到解析解,则返回其积分形式

省略自变量和指定自变量的例子:

>> dsolve('Dy = 2*x' , 'x')    %自变量x
ans =
x^2 + C1
 
>> dsolve('Dy = 2*x')    %默认自变量t
ans =
C2 + 2*t*x

解析解例题:

例1:求微分方程 dy/dx+2xy=xe^{-x^{2}}通解

>> y=dsolve('Dy+2*x*y = x*exp(-x^2)' , 'x')
y =
C3*exp(-x^2) + (x^2*exp(-x^2))/2

例2:求下面微分方程的特解

>> y=dsolve('D2y+4*Dy+29*y=0' , 'y(0)=0 ,Dy(0)=15' , 'x')
y =
3*sin(5*x)*exp(-2*x)

例3:求微分方程组的通解

>> [x , y , z]=dsolve('Dx =2*x-3*y+3*z' , 'Dy=4*x-5*y+3*z' , 'Dz=4*x-4*y+2*z' , 't')

x =
C6*exp(2*t) + C7*exp(-t)

y =
C6*exp(2*t) + C7*exp(-t) + C8*exp(-2*t)

z =
C6*exp(2*t) + C8*exp(-2*t)

二、数值解

只有很少一部分微分方程 (组) 能求出解析解;

大部分微分方程(组)只能利用 数值方法 近似求解

1、solver函数:

  • [T , Y] = solver(odefun , tspan , y0)
    • odefun——待求解的微分方程,可用matlab函数表示(m文件名)
    • tspan——求解区间
    • y0——初值条件(即自变量取最小值时(不一定为0)的函数值)
    • T——(向量)返回分割点的值,作自变量
    • Y——(向量)返回函数在分割点上的函数值
  • MATLAB在数值求解时 自动对求解区间进行分割
  • solver 为MATLAB的 ODE求解器(求解器中可以调用各种数值求解法:ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb)
求解器特点说明
ode454,5阶 R-K 法首选方法
ode232,3阶 R-K 法精度较低
ode113Adams法,高低精度均可到 10^-3 ~10^-6计算时间比ode45短
ode23t梯形算法(改进欧拉法)
ode15s多步法:精度中等ode45失效时,可用
ode23s单步法:低精度当精度较低时,计算时ode15s短
ode23tb梯形算法:低精度当精度较低时,计算时比ode15s短

——————————————————————————————————

2、odenfun(待求解的微分方程)的构建:

(1)odenfun为 显示常微分方程(类似于显函数,可分离变量)

  • 可以用 inline函数 定义:直接将函数表达式内嵌在命令行里
    变量名=inline('函数表达式' , '变量名1' , '变量名2' , ...)
  • 也可以在 函数文件 中定义,然后通过函数句柄调用

例1:一阶常微分非线性方程数值解,求解范围 [0,0.5]

>> fun=inline('-2*y+2*x^2+2*x','x','y');	%定义常微分方程
>> [x,y] = ode23(fun , [0,0.5] , 1)		%数值求解法使用ode23,自变量范围[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
    
    
    
    
>> [x,y] = ode23(fun , [0:0.1:0.5] , 1)		%自定义步长为0.1
x =
         0
    0.1000
    0.2000
    0.3000
    0.4000
    0.5000

y =
    1.0000
    0.8287
    0.7103
    0.6388
    0.6093
    0.6179

——————————————————————————————————————

(2)求 高阶常微分方程,需将其化为 一阶常微分方程组 ,用函数文件定义

函数文件格式 function dy = vdp1000(x , y)————其中x是向量,y是

数值解例题:

例:求解:

解: 令 

则原式子化成:

vdp1000.m

function dy = vdp1000(t , y)	%编写函数vdp1000,输入的t是自变量(不可省略),在后面的数值计算时,将自变量范围赋给t;y是一个2列的行向量,因为我们后面赋初值是有2个,且输出值也有2个
dy = zeros(2,1);		%建立一个2行的列向量
dy(1) = y(2);			%按照上边的式子赋值
dy(2) = 1000*(1-y(1)^2)*y(2)-y(1);		%按照上边的式子赋值

命令窗口

>> [T , Y]=ode15s('vdp1000' , [0 , 3000] , [2 0]);		%调用方法,自变量t的范围[0 , 3000],将初值 2,0 赋给 y(1),y(2) ,T——n个元素的列向量,Y——n行2列的矩阵,Y(:1)是式子中的y1,Y(:2)是式子中的y2
>> plot(T , Y(:,1))		%显示Y(:1),即式子中的y1,即原始式子的因变量x

例2:手写 欧拉法 求常微分方程,x 范围 [0,0.5] , 步长 h=0.1

解:根据欧拉公式: yn+1=yn+hf(xn,yn) -------------------(注:这里的 f(x,y)=y′ )

代入 y′ , 可以得到 yn+1=0.9yn+0.1xn+0.1

Euler.m

function y=Euler()

x=[0:0.1:0.5];		%定义x
y=zeros(size(x));	%创建y
h=0.1;				%步长h=0.1
len = length(x);	%向量x的长度
y(1)=1;				%y初始化
for i=2:len			%从y2开始计算
    y(i)=0.9*y(i-1)+h*x(i-1) + 0.1;
end

%画图
figure
hold on 
grid on
plot(x , y ,'r')		%红色线为数值解
plot(x , x+exp(-x),'b')	%蓝色线为解析解
xlabel('x')
ylabel('y')
title('Euler')

命令窗口

%函数调用
>> y=Euler()
y =
    1.0000    1.0000    1.0100    1.0290    1.0561    1.0905

例3:追击问题

设位于坐标原点的甲舰向位于 x 轴上 A(1,0) 处的乙舰发射导弹,导弹头始终对准乙舰。

已知:乙舰以最大的速度 v0 沿平行于 y 轴的直线行驶,导弹的速度是 5v0 。

要求:模拟导弹的运行轨迹,并求击中点的坐标

解:

建立坐标系

设 导弹的飞行曲线方程 为 y=y(x)

在时刻 t 位于 P(x,y)

eql.m

function dy=eql(x,y)	%一阶微分方程组
dy=zeros(2,1);		%构造2行的列向量
dy(1)=y(2);				%dy(1)
dy(2)=1/5*sqrt(1+y(2)^2)/(1-x);		dy(2)

命令串口

>> [x , y]=ode15s('eql',[0:0.001:1],[0 0]);		%数值求解微分方程组
>> plot(x,y(:,1),'r.')					%画导弹轨迹,dy(1),即原式子的y
hold on
>> y2=0:0.001:2;								
>> plot(1,y2,'b*')					%画乙舰轨迹
axis( [0 , 1 , 0 , 0.25 ] )				%坐标轴范围
>> legend('导弹轨迹' , '乙舰轨迹')			%图例

(二)常微分方程的MATLAB求解——手写龙格-库塔法(R-K法)模板、一般步骤 

概要:

  • 龙格-库塔(R-K)法的写法:就是不断调用微分方程组,迭代计算出对于K1,K2,...,最后再叠加。需要注意的是高阶微分方程,其原函数的导数也是通过迭代计算得到的
  • 本人归纳了其套用 R-K 法的一般套路:3个函数、3个步骤——这也是MATLAB自带的求解方法的步骤
    • 三个函数:
      • Fun函数——用于存放一阶微分方程组
      • RK函数——用于使用几阶的R-K法求数值解,上边我只写了 4阶R-K法
      • 赋初值函数——只是单纯的写 x的范围,步长h,矩阵y的阶数,原函数的各个初值;以及调用 RK函数
    • 三个步骤:
      • 赋初值:写 x的范围,步长h,矩阵y的阶数,原函数的各个初值
      • 将高阶微分方程 拆分成 一阶微分方程组
      • 修改 Fun 函数:注意——向量dy的长度,和高阶微分方程的阶次有关
  • 本人也提供了 手写R-K法的模板,供大家使用参考

例一:一阶微分方程

使用MATLAB自带的 ode45函数 和 手写4阶R-K法 求解下列方程

1.ode45函数解法:

微分方程组

function dy=Fun(x,y)
dy=zeros(size(y));

dy(1) = sin(x)+y(1);					   %dy(1)表示以y的一阶导为f(x,y)

调用ode45函数

function y = MATLAB_RK()
[X ,Y]=ode45('Fun' , [0 :0.1 :20] , [1]);		% y的初值1

%画图
hold on
grid on
plot(X , Y(:,1))		%输出 Y的第一列,即原函数关于x的曲线
xlabel('x')
ylabel('y')

2.手写4阶R-K解法:

微分方程组

function dy=Fun(x,y)
dy=zeros(size(y));

dy(1) = sin(x)+y(1);					   %dy(1)表示以y的一阶导为f(x,y)

调用手写R-K函数

function y=PlotAll()
% 准备阶段
x=[0 : 0.1 :20];			% x——行向量
h=0.1;					% 步长
y=zeros(length(x) , 1);		        % y——矩阵,行数是x的元素个数,列数是几阶微分方程
					% 后面的1表示一阶微分方程,如果是2则表示二阶微分方程,以此类推
y(1,:)=[1];      			% y,第一行赋初值,几阶微分方程就要赋几个初值

% 函数调用
Y=FourRK(x,h,y);			% 调用4阶R-K法,输入参数x,h,y


%画图
hold on
grid on
plot(x , Y(:,1))		%输出 Y的第一列,即原函数关于x的曲线
xlabel('x')
ylabel('y')

4阶R-K函数

function y=FourRK(x,h,y)

len = length(x);	

for i=2:len		%循环:直到求完最后一个x取值

K1 = Fun(x(i-1),y(i-1,:));  		% K1

K2 = Fun(x(i-1)+1/2*h , y(i-1,:)+1/2*h.*K1);		% K2

K3 = Fun(x(i-1)+1/2*h , y(i-1,:)+1/2*h.*K2);		% K3

K4 = Fun(x(i-1)+h , y(i-1,:)+h.*K3);			% K4


y(i,:) = y(i-1,:)+h/6.*(K1 + 2*K2 + 2*K3 + K4);           % 更新y的第i行,y(i,1)是原函数的y,y(i,2)是原函数的一阶导,以此类推
														  % 因为这里是1阶微分方程,所以只有y(i,1),无y(i,2)
end

例二:二阶微分方程

使用MATLAB自带的 ode45函数 和 手写4阶R-K法 求解下列方程

解:将高阶微分方程 化成 一阶微分方程组

我们下面的 Fun函数 就是依据这个一阶微分方程组写的

1.ode45函数解法:

微分方程组

function dy=Fun(x,y)			  % y(1)表示原函数y,y(2)表示原函数y的一阶导
dy=zeros(size(y));
dy(1)=y(2);                               % dy(1)表示以y的一阶导为f(x,y)
dy(2)=-1/2*(y(1)+(y(1)*y(1)-1)*y(2));     % dy(2)表示以y的二阶导为f(x,y)

调用ode45函数

function y = MATLAB_RK()
[X ,Y]=ode45('Fun' , [0 :0.1 :20] , [0.25 , 0.0]);		% y的初值0.25 ,y一阶导的初值 0.0

%画图
hold on
grid on
plot(X , Y(:,1))		%输出 Y的第一列,即原函数关于x的曲线
xlabel('x')
ylabel('y')

2.手写4阶R-K解法:

微分方程组

function dy=Fun(x,y)
dy=zeros(size(y));
dy(1)=y(2);                               %dy(1)表示以y的一阶导为f(x,y)
dy(2)=-1/2*(y(1)+(y(1)*y(1)-1)*y(2));     %dy(2)表示以y的二阶导为f(x,y)

调用手写R-K函数

function y=PlotAll()
% 准备阶段
x=[0 : 0.1 :20];		% x——行向量
h=0.1;				% 步长
y=zeros(length(x) , 2);	        % y——矩阵,行数是x的元素个数,列数是几阶微分方程
				% 后面的1表示一阶微分方程,如果是2则表示二阶微分方程,以此类推。这里是二阶微分方程,所以是2
y(1,:)=[0.25 , 0.0];            % y,第一行赋初值,几阶微分方程就要赋几个初值。
				% y的初值0.25 ,y一阶导的初值 0.0

% 函数调用
Y=FourRK(x,h,y);		% 调用4阶R-K法,输入参数x,h,y

%画图
hold on
grid on
plot(x , Y(:,1))		%输出 Y的第一列,即原函数关于x的曲线
xlabel('x')
ylabel('y')

4阶R-K函数

function y=FourRK(x,h,y)

len = length(x);

for i=2:len			%循环:直到求完最后一个x取值

K1 = Fun(x(i-1),y(i-1,:));  		% K1


K2 = Fun(x(i-1)+1/2*h , y(i-1,:)+1/2*h.*K1);		% K2


K3 = Fun(x(i-1)+1/2*h , y(i-1,:)+1/2*h.*K2);		% K3


K4 = Fun(x(i-1)+h , y(i-1,:)+h.*K3);		% K4


y(i,:) = y(i-1,:)+h/6.*(K1 + 2*K2 + 2*K3 + K4);           % 更新y的第i行,y(i,1)是原函数的y,y(i,2)是原函数的一阶导,以此类推
														  % 因为这里是2阶微分方程,所以有y(i,1),有y(i,2)
end





%本质上,将高阶微分化为一阶微分方程组,然后按照各自当作f(x,y),分别进行迭代,最后更新

例二 手写R-K法 答疑:

一阶微分方程可能还好,二阶微分方程可能就没那么好理解了,为此对于例二我换了种 不通用但更直观 的写法,帮助大家理解

首先得理解为什么要将 二阶微分方程 拆成 一阶微分方程组 ,因为这样我们就可以根据 设 f(x,y)=y 的一阶导数,进行迭代求解。————如此题,我们拆成了2个一阶微分方程,所以就需要就同时进行2个R-K迭代

  • 这里我将 y 和 y' 拆开成2个行向量——————上边是 将2个合并到一个矩阵 y 中, y(1)为原函数y,y(2)为原函数y的一阶导
  • 并且我们需要分别定义函数,来迭代计算 y 和 y' ——————上边是 将2个合并到一个矩阵 y 中.一起计算了

Fun1————迭代求 y

function f1=Fun1(x,y,y1)

f1=y1;      %f2表示以y的一阶导为f(x,y)——f(x,y)=输入的y1

Fun2————迭代求 y‘

function f2=Fun2(x,y,y1)

f2=-1/2*(y+(y*y-1)*y1);     %f2表示以y的二阶导为f(x,y)

4阶R-K函数

function y=FourRK2()

x=[0 : 0.1 :20];
y=zeros(size(x));
y1=zeros(size(x));
h=0.1;
len = length(x);
y(1)=0.25;      % 原函数y初值
y1(1)=0.0;      % 原函数y的一阶导初值

for i=2:len
%Kmn 表示第n次迭代,以m次导数作为f(x,y)

K11 = Fun1(x(i-1),y(i-1),y1(i-1));  
K21 = Fun2(x(i-1),y(i-1),y1(i-1));

K12 = Fun1(x(i-1)+1/2*h , y(i-1)+1/2*h*K11 , y1(i-1)+1/2*h*K21);
K22 = Fun2(x(i-1)+1/2*h , y(i-1)+1/2*h*K11 , y1(i-1)+1/2*h*K21);

K13 = Fun1(x(i-1)+1/2*h , y(i-1)+1/2*h*K12 , y1(i-1)+1/2*h*K22);
K23 = Fun2(x(i-1)+1/2*h , y(i-1)+1/2*h*K12 , y1(i-1)+1/2*h*K22);

K14 = Fun1(x(i-1)+h , y(i-1)+h*K13 , y1(i-1)+h*K23);
K24 = Fun2(x(i-1)+h , y(i-1)+h*K13 , y1(i-1)+h*K23);

y(i) = y(i-1)+h/6*(K11 + 2*K12 + 2*K13 + K14);           % 更新下一个原函数y
y1(i) = y1(i-1)+h/6*(K21 + 2*K22 + 2*K23 + K24);         % 更新下一个原函数y的一阶导
end


%画图
hold on
grid on
plot(x , y)		
xlabel('x')
ylabel('y')


%本质上,将高阶微分化为一阶微分方程组,然后按照各自当作f(x,y),分别进行迭代,最后更新

总结:手写微分方程的 模板 和 一般步骤

一、原理讲解:

1.模板:

  • Fun函数——用于存放一阶微分方程组
  • RK函数——用于使用几阶的R-K法求数值解,上边我只写了 4阶R-K法
  • 赋初值函数——只是单纯的写 x的范围,步长h,矩阵y的阶数,原函数的各个初值;以及调用 RK函数

2.一般步骤:

其实这里 手写的R-K法 原理和步骤就是 MATLAB自带的求解R-K方法

  • 赋初值:写 x的范围,步长h,矩阵y的阶数,原函数的各个初值
  • 将高阶微分方程 拆分成 一阶微分方程组
  • 修改 Fun 函数:注意——向量dy的长度,和高阶微分方程的阶次有关

二、具体模板:

1.微分方程模板和准备函数模板

Fun函数————需要自己修改

function dy=Fun(x,y)
dy=zeros(size(y));
dy(1)=    ;                               %dy(1)表示以y的一阶导为f(x,y)
dy(2)=    ;     			  %dy(2)表示以y的二阶导为f(x,y)
dy(3)=    ;     			  %dy(3)表示以y的三阶导为f(x,y)
...

准备初值函数————需要自己修改值,需要修改的用【】括起了

function Y=PlotAll()
% 准备阶段
x=[0 : 0.1 :20];		% 【x——行向量】
h=0.1;				% 【步长】
y=zeros(length(x) , 2);	        % 【y——矩阵,行数是x的元素个数,列数是几阶微分方程】
				% 后面的1表示一阶微分方程,如果是2则表示二阶微分方程,以此类推。这里是二阶微分方程,所以是2
y(1,:)=[0.25 , 0.0];            % 【y,第一行赋初值,几阶微分方程就要赋几个初值。】
				% 【y的初值0.25 ,y一阶导的初值 0.0】

% 函数调用
Y=FourRK(x,h,y);		% 调用4阶R-K法,输入参数x,h,y

%画图
hold on
grid on
plot(x , Y(:,1))		%输出 Y的第一列,即原函数关于x的曲线
xlabel('x')
ylabel('y')

2.欧拉法、改进欧拉法、4阶R-K法模板

1阶R-K函数(欧拉法)

function y=OneRK(x,h,y)

len = length(x);

for i=2:len			%循环:直到求完最后一个x取值

K1 = Fun(x(i-1),y(i-1,:));

y(i,:)=y(i-1,:)+h.*K1;		% 更新y的第i行,y(i,1)是原函数的y,y(i,2)是原函数的一阶导,以此类推
	
end





%本质上,将高阶微分化为一阶微分方程组,然后按照各自当作f(x,y),分别进行迭代,最后更新

2阶R-K函数(改进欧拉法)

function y=TwoRK(x,h,y)

len = length(x);

for i=2:len			%循环:直到求完最后一个x取值

K1 = Fun(x(i-1) , y(i-1,:));
K2 = Fun(x(i-1)+h , y(i-1,:)+h.*K1);

y(i,:)=y(i-1,:)+h/2.*(K1+K2);		% 更新y的第i行,y(i,1)是原函数的y,y(i,2)是原函数的一阶导,以此类推
end





%本质上,将高阶微分化为一阶微分方程组,然后按照各自当作f(x,y),分别进行迭代,最后更新

4阶R-K函数

function y=FourRK(x,h,y)

len = length(x);

for i=2:len			%循环:直到求完最后一个x取值

K1 = Fun(x(i-1),y(i-1,:));  		% K1


K2 = Fun(x(i-1)+1/2*h , y(i-1,:)+1/2*h.*K1);		% K2


K3 = Fun(x(i-1)+1/2*h , y(i-1,:)+1/2*h.*K2);		% K3


K4 = Fun(x(i-1)+h , y(i-1,:)+h.*K3);		% K4


y(i,:) = y(i-1,:)+h/6.*(K1 + 2*K2 + 2*K3 + K4);           % 更新y的第i行,y(i,1)是原函数的y,y(i,2)是原函数的一阶导,以此类推
end



%本质上,将高阶微分化为一阶微分方程组,然后按照各自当作f(x,y),分别进行迭代,最后更新

补充:关于解为lambertw的微分方程数值解

例,求解下面的微分方程

可以发现这个初值会让 y′(0) 的值有无数个。如果我们写 dy=xy/(1−y) 就会出现分母为0的情况,一开始我就是这么定义 Fun 函数的,发现结果为 NaN。

正解:

  • 其实这个函数与 y′ 没有关系,不管 y′(0) 取多少,函数图像总是先为0,后面快速趋向正无穷
  • 因此我们就默认取 dy=0 (当然大家也可以取其他值),然后不将 dy 到一边,直接用原表达式计算

Fun

function dy=Fun(x,y)

dy=-200;
dy=(x*y)+dy*y;

Logo

开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!

更多推荐