用MATLAB求解微分方程

2024-06-26 18:47| 来源: 网络整理| 查看: 265

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

（一）常微分方程的MATLAB求解概要：常微分方程的MATLAB求解分为解析解、数值解解析解(只有少数微分方程组有解析解)：dsolve函数数值解：solver函数，MATLAB内部自带了许多求解器（各种数值求解的方法，如 R-K法）常用函数文件构造微分方程组数值求解时：高阶微分方程化成常微分方程组一、解析解

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;

【本文地址】

公司简介

联系我们