MATLAB-求解常微分方程

最新推荐文章于 2024-04-23 17:09:42 发布

hugehugg

最新推荐文章于 2024-04-23 17:09:42 发布

阅读量1.5k

点赞数

文章标签： matlab

本文链接：https://blog.csdn.net/hugehugg/article/details/132083812

版权

用matlab方法求解常微分方程

欧拉法

对于一阶常微分方程

$\begin{cases}y'=f(x,y)\\y(x_0)=y_0\end{cases}$

对应的改进的欧拉公式为

$\begin{cases}y_p=y_n+hf(x_n,y_n)\\y_q=y_n+hf(x_n+h,y_p)\\y_{n+1}=\dfrac{1}{2}(y_p+y_q)\end{cases}$

对应的matlab代码为

function [x,y]=euler(fun,x0,xfinal,y0,n)
%检查是否有足够的输入参数值传入函数。
%如果输入参数的数量小于 5，说明用户没有提供 n 的值，
%那么将默认将 n 设置为 50。
if nargin<5
    n=50;
end

h=(xfinal-x0)/n;%步长
x(1)=x0;y(1)=y0;%初始条件

%使用函数句柄来表示和传递函数
%result = feval(function_handle, arg1, arg2, ...)
%function_handle 是一个函数句柄，arg1、arg2 等是要传递给该函数的参数
for i=1:n
    x(i+1)=x(i)+h;
    y1=y(i)+h*feval(fun,x(i),y(i));
    y2=y(i)+h*feval(fun,x(i+1),y(i));
    y(i+1)=(y1+y2)/2;
end

例子

$y'=-2y+2x^{2}+2x,\left(0\leq x\leq0.5\right),y(0)=1$

对应代码

function f=fun(x,y)
f=-2*y+2*x^2+2*x;
end

在matlab命令行窗口输入

[x,y]=euler('fun',0.5,1,10)

得到数值解

函数功能介绍及例子

dsolve函数

用于求解常微分方程组的精确解。若没有初始条件或边界条件，则求出通解；若有，则求出特解。

函数格式

Y=dsolve('eq1,eq2,...','cond1,cond2,...','name');
%eq1:微分方程
%cond1：初始条件或边界条件
%name：表示变量。若没有变量，默认为t

例题

$\frac{dy}{dx}=3x^2.$

例题（加上初始条件）

$\frac{dy}{dx}=3x^2.$ $y(0)=2$

>> dsolve('Dy=3*x^2','x')
 
ans =
 
x^3 + C2

>> dsolve('Dy=3*x^2','y(0)=2','x')
 
ans =
 
x^3 + 2

ode函数（ode23,ode45,ode113使用方法）

符号介绍

D：微分符号；D2：二阶微分；D3：三阶微分

上文中的dsolve函数存在大量常微分方程。从理论上讲，解是存在的但无法求出解析解（就i是一些严格的公式）。因此，但现在要寻求数值解（利用有限元或插值等方法求得的近似解）。

《数值分析》

例如求解微分方程 $\frac{dy}{dx}=3x^2.$ 可以转化为 $\frac{y(n)-y(n-1)}{dt}=3t^{2}$ .即 $y(n)=y(n-1)+3t^{2}*dt$ .就可以用迭代的方式求解y。dt看作是步长。

ode是matlab专门用于求解微分方程的功能函数。该求解器有变步长和定步长两种类型。

ode45	4、5阶R-K方程	首选，没用的话换ode15s
ode23	2、3阶R-K方程	精度较差
ode113	多部算法	计算时间短

下面以ode45为例介绍用法

函数格式

[T,Y]=ode45('odefun',tspan,y0)

[T,Y]=ode45('odefun',tspan,y0,options)

[T,Y,TE,YE,IE]=ode45('odefun',tspan,y0)

sol=ode45('odefun',[t0 tf],y0...)

其中：

ode是函数句柄，函数文件名，匿名函数句柄，内联函数名

tspan是求解区间[t0 tf]

y0是初始值向量

T是返回列向量的时间点

Y是对应T的求解列向量

options是求解参数设置

TE 事件发生事件

YE 事件发生时的答案

IE事件函数消失时的指针i

例题：

$y^{\prime}=-2y+2x^{2}+2x,\left(0\leq x\leq0.5\right),y(0)=1$

function f=doty(x,y)
f=-2*y+2*x^2+2*x;

%%命令行窗口输入
[x,y]=ode45('doty',0,0.5,1)

微分方程标准化

利用ode45求解高阶微风方程时，需要做变量替换。

微分方程为 $F(y,y',y'',...,y^{(n)},t)=0$

初始条件 $y(t_{0})=c_{0},y^{\prime}(t_{0})=c_{1},...,y^{(n-1)}(t_{0})=c_{n-1}$

首先做变量替换 $x(1)=y,x(2)=y',x(3)=y'',...,x(n)=y^{(n-1)}.$

原微分方程可以转换为下面的微分方程组格式：

$\begin{cases}d\text{x}(1)=x(2)\\d\text{x}(2)=x(3)\\......\\d\text{x}(n-1)=x(n)=f(t,x(n-1),...x(2),x(1))\end{cases}$

下面就可以利用转化好的微分方程组来编写odefun函数

例题：求解常微分 $y^{\prime\prime}=-ty+e^{t}y^{\prime}+3\sin2t$ 。

先编写odefun函数：

令 $x(1)=y,x(2)=y\prime$

则微分方程可以转换为：

$\begin{aligned}dx(1)&=y'=x(2)\\dx(2)&=y''=-t^*x(1)+\exp(t)*x(2)+3*\sin(2*t)\end{aligned}.$

function dx=odefun(t,x)
dx=zeros(2,1); %初始化dx
dx(1)=x(2);
dx(2)=-t*x(1)+exp(t)*x(2)*sin(2*t);
end

例题：

在区间[3,9,4]内求解常微分方程 $y"=-ty+e^{t}y'+3\sin2t,$ $y(t=3)=8,y'(t=3)=2$ .

在上面的基础上编写主函数，加上求解区间和边值条件。

注意ode的运行结果以列向量形式给出。在本例中，x的第一列为y,第二列为y‘，若遇到变量不是列向量，可以考虑用reshape函数做矩阵变换。

则，plot(t,x(:,1))画出来是x的第一列数据，即y。plot(t,x(:,2))画出来是x的第二列数据，即y’。

tspan=[3.9,4];%求解区间
y0=[8 2];%初值
[t,x]=ode45('odefun',tspan,y0)

%%画图
plot(t,x(:,1),'-o',t,x(:,2),'-*')
xlabel('t')
ylabel('y')

解一般形式常微分方程的通解

无初值无边界的常微分方程的解就是方程的通解

dsolve('diff_equ')  %diff_equ为常微分方程组，默认t为自变量
dsolve('diff_equ'，‘var’)  %diff_equ为常微分方程组，var为自定义变量

例题： $x^2+y+(x-2y)y^{\prime}=0$

syms x y
diff_equ='x^2+y+(x-2*y)*Dy=0';
dsolve(diff_equ,'x')

求解线性常微分方程组

一阶齐次线性微分方程组

$X'=AX,X=\begin{bmatrix}x_1\\\vdots\\x_n\end{bmatrix},A=\begin{bmatrix}a_{11}&\cdots&a_{1n}\\\vdots&\vdots&\vdots\\a_{n1}&\cdots&a_{nn}\end{bmatrix}$ , $X(t_0)=X_0$ 的解为 $X(t)=e^{A(t-t_0)}X_0$ 。

例题： $X'=\begin{bmatrix}2&1&3\\0&2&-1\\0&0&2\end{bmatrix}X,X(0)=\begin{bmatrix}1\\2\\1\end{bmatrix}$

syms t
a=[2,1,3;0,2,-1;0,0,2];
x0=[1;2;1];
x=expm(a*t)*x0

一阶非齐次

$X'=AX+f(t),X(t_0)=X_0$

解为 $X(t)=e^{A(t-t_0)}X_0+\int_{t_0}^te^{A(t-s)}f(s)ds.$

例题： $X'=\begin{bmatrix}1&0&0\\2&1&-2\\3&2&1\end{bmatrix}X+\begin{bmatrix}0\\0\\e'\cos2t\end{bmatrix},X(0)=\begin{bmatrix}0\\1\\1\end{bmatrix}$

clc,clear
syms t s
a=[1,0,0;2,1,-2;3,2,1];ft=[0;0;exp(t)*cos(2*t)];
x0=[0;1;1];
x=expm(a*t)*x0+int(expm(a*(t-s))*subs(ft,s),s,0,t);

hugehugg

关注

0
点赞
踩
4

收藏

觉得还不错? 一键收藏
0
评论
MATLAB-求解常微分方程

用于求解常微分方程组的精确解。若没有初始条件或边界条件，则求出通解；若有，则求出特解。函数格式%eq1:微分方程%cond1：初始条件或边界条件%name：表示变量。若没有变量，默认为t例题例题（加上初始条件）ans =x^3 + C2ans =x^3 + 2D：微分符号；D2：二阶微分；D3：三阶微分上文中的dsolve函数存在大量常微分方程。从理论上讲，解是存在的但无法求出解析解（就i是一些严格的公式）。因此，但现在要寻求数值解（利用有限元或插值等方法求得的近似解）。
复制链接

扫一扫