大家好,又见面了,我是你们的朋友全栈君。
上篇博客介绍了Matlab求解常微分方程组解析解的方法:博客地址 微分方程组复杂时,无法求出解析解时,就需要求其数值解,这里来介绍。 以下内容按照Matlab官方文档提供的方程来展开(提议多看官方文档)
一般形式:[t,y] = ode45(odefun,tspan,y0)
其中 tspan = [t0 tf]
功能介绍:求微分方程组 y′=f(t,y) 从 t0 到 tf 的积分,初始条件为 y0。解数组 y 中的每一行都与列向量 t 中返回的值相对应。
方程: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
中间是空格,不能用逗号
结果:
方程:范德波尔方程 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’’,然后把两个方程写成列向量形式就行了 结果:
这里就和求解二阶方程类似的,只不过不需要降阶,仍旧需要一个函数来定义方程组。我们这里不用官方文档的例子,用同学的循坏摆问题来进行演示。 方程:
给定的初值(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即可。
结果:
讲到这里,大部分我们用到的微分方程形式都可以求解了,Matlab还支持带有时变项和额外参数的微分方程求解,这里不再赘述,大家可以自行参阅官方文档。
发布者:全栈程序员栈长,转载请注明出处:https://javaforall.cn/126077.html原文链接:https://javaforall.cn