Matlab求常微分方程组的数值解

上篇博客介绍了Matlab求解常微分方程组解析解的方法:博客地址
微分方程组复杂时,无法求出解析解时,就需要求其数值解,这里来介绍。
以下内容按照Matlab官方文档提供的方程来展开(提议多看官方文档)

介绍一下核心函数ode45()

一般形式:[t,y] = ode45(odefun,tspan,y0) 其中 tspan = [t0 tf]
功能介绍:求微分方程组 y′=f(t,y) 从 t0 到 tf 的积分,初始条件为 y0。解数组 y 中的每一行都与列向量 t 中返回的值相对应。

1. 一阶微分方程求解(简单调用即可)

方程:y’=2*t
代码

tspan=[1 6]; %定义自变量x的取值空间为1-6
y0=0;%定义因变量的初值,当x=1(x取值空间的第一个数)时,y0=0
[t,y]=ode45(@(t,y) 2*t,tspan,y0); %定义函数y'=2*t,使用ode45求解
plot(t,y,'-o'); %绘制求得的数值曲线

说明:简单的odefun参数就是这个形式,@(x,y) fun 中间是空格,不能用逗号
结果
在这里插入图片描述

2. 二阶微分方程求解(引入函数文件)

方程:范德波尔方程 y1’’-u(1-y1²)*y1’+y1=0;(这里设u=1)
代码
定义输入的方程,以函数形式定义

function dydt=odefun(t,y)
%二阶方程为y1''-(1-y1²)*y1'+y1=0;
%降阶为两个方程:y1'=y2; 
%               y2'=(1-y1²)*y2-y1;
%t虽然没有使用,但必须要作为参数写入
dydt=[y(2);(1-y(1)^2)*y(2)-y(1)];

求解作图

tspan=[0 20]; %定义自变量x的取值空间为0-20
y0=[2;0];%定义因变量的初值,当x=0时,y1=2,y2=y1'=0;
[t,y]=ode45(@odefun,tspan,y0); %使用ode45求解
%%下面为作图过程,不解释
plot(t,y(:,1),'-o',t,y(:,2),'-o')
title('Solution of van der Pol Equation (\mu = 1) with ODE45');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')

说明:先进行降阶,使y2=y1’,这时候y2’=y1’’,然后把两个方程写成列向量形式就行了
结果
在这里插入图片描述

3. 求解微分方程组(和2类似)

这里就和求解二阶方程类似的,只不过不需要降阶,仍旧需要一个函数来定义方程组。我们这里不用官方文档的例子,用同学的循坏摆问题来进行演示。
方程在这里插入图片描述
给定的初值(w接近0,但实际上不能设置为0):
在这里插入图片描述
代码
定义输入的方程

function dRvw=func(t,Rvw)
%% 函数功能:为ode45提供微分方程
%输入:t:时间序列,就是θ;Rvw:因变量,Rvw(1)代表R,Rvw(2)代表v,Rvw(3)代表w
%输出:dRvw:因变量的一阶微分,dRvw(1)代表dR,dRvw(2)代表dv,dRvw(3)代表dw
%% 初始化因变量的一阶微分,3×1的向量
dRvw=zeros(3,1);
%% 参数初始化
r=0.01;u=0.1;g=9.8;M=10;m=1;
%% 输入微分方程式
dRvw(1)=-Rvw(2)/Rvw(3)-r;
dRvw(2)=(M*g-(m*Rvw(3)^2*Rvw(1)+m*g*sin(t)) *exp(u*(t+pi/2)))/(Rvw(3)*(M+m));
dRvw(3)=Rvw(3)*r/Rvw(1)+2*Rvw(2)/Rvw(1)+g*cos(t)/(Rvw(3)*Rvw(1));

求解作图

%% 参数初始化 start_Theta是θ的初始值 end_Theta是θ的结束值
%R是半径初值;v是线速度初值;w是角速度初值
start_Theta=0;end_Theta=2*pi;R=1;v=0;w=1e-5;
%% 使用ode45方法计算微分方程组func的数值解
%func是带有方程组的函数 
%[start_Theta end_Theta]是自变量范围
%[R;v;w]是方程初值
%T是自变量的数组,Rvw是对应的因变量的数值。Rvw(:,1)=R;Rvw(:,2)=v;Rvw(:,3)=w;
[T,Rvw]=ode45(@func,[start_Theta end_Theta],[R;v;w]);
%% 组合自变量和因变量。TRvw(:,1)=θ;TRvw(:,2)=R;TRvw(:,3)=v;Trvw(:,4);
TRvw=[T,Rvw];
%% 绘制自变量-因变量图像.figure1是R-θ,figure2是v-θ,figure3是w-θ
plot(T,Rvw(:,1),'-o',T,Rvw(:,2),'-o',T,Rvw(:,3),'-o')
title('自变量-因变量图像');
xlabel('Angle θ');
ylabel('Solution');
legend('R','v','w')

说明:注释的应该是比较清楚的,把三个方程写成列向量的形式就行

PS:有些人和我说不能运行,然后我看了他们出错的原因,有点儿哭笑不得。出错的基本上都是运行上面的dRvw=func(t,Rvw)这个函数的。说明一下,这是有参数的函数,不给参数不能直接运行的。下面的求解作图脚本才是需要运行的哈,它调用了函数,才得到的结果。如果大家还发现什么问题,欢迎私戳或评论。

PS+ 有了PS之后,还是很多人问我参数的问题,我在这里直接把文件给大家:cupt.zip
提取码:6k8n。
解压后,有两个文件,大家直接运行cupt.m即可。

结果
在这里插入图片描述

4. 更多形式

讲到这里,大部分我们用到的微分方程形式都可以求解了,Matlab还支持带有时变项和额外参数的微分方程求解,这里不再赘述,大家可以自行参阅官方文档。

已标记关键词 清除标记
相关推荐
©️2020 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页