MATLAB绘制微分方程的相图/方向场/向量场

对于微分方程

$$y^{\prime }=1-y^2$$

画出他的相图，其中 $x\in [-3,3]$ ,$y\in [-3.3]$,步长为 $h=0.2$,箭头长为 $c=0.01$

则MATLAB代码为

clc,clear,close all
x_0=-3:0.2:3;
y_0=-3:0.2:3;
[x,y]=meshgrid(x_0,y_0);
d=sqrt(1+(1-y.^2).^2);
u=1./d;
v=1*(1-y.^2)./d;
quiver(x,y,u,v);
xlim([-3,3])
ylim([-3,3])
%		下面几行代码用来绘制初值为 y(0)=1.6 的数值解图像
%		[X,Y]=ode45(@(x,y) 1-y.^2,[-3,3],[0,1.6]);
%		hold on
%		plot(X,Y(:,2),'-r',"LineWidth",2)
%		hold off
%%%%%%%%%% 以下代码保存成 ode1.m 文件 %%%%%%%%%%
%		function dy= ode1(~,y)
%    		dy = 1:1-y(2).^2;
%    	end

---

2021年2月9日20:35:54

---

有朋友问 u v 是怎么算的，现做如下解释。

我 quiver 函数用的不多，存在用错的可能。如果我的解释有误，还望指教。

首先我们要知道画向量场的原理是什么。画向量场其实就是画很多箭头。那么为了确定箭头的形式，要确定箭头的起点，箭头的方向，箭头的大小。

箭头的起点，这里采用是在 $[-3,3]\times [-3,3]$ 上等距取点，也就是第二三四行代码。

那么箭头的方向和大小如何计算呢？首先要知道箭头的方向是什么含义。箭头的方向表示曲线在该点处的切线方向，而一点处的切线方向，可以用 $\Delta x,\Delta y$ 来表示。

而我们要算的 u v，就是 $\Delta x,\Delta y$。

由微积分的知识我们可以知道。

$$\frac{\Delta y}{\Delta x}=y^{\prime }\text{\hspace{1em}(1)}$$

可以想象到在该切线方向上，有任意多满足 $(1)$ 的 $(\Delta x,\Delta y)$ ，但我们这里考虑单位向量，即：

$$\Delta x^2+\Delta y^2=1\text{\hspace{1em}(2)}$$

至于我为什么要考虑单位向量呢？原因有二，一个是总要有一个标准来确定 $\Delta x,\Delta y$，选啥向量都一样，还不如选个单位向量感觉比较酷。另一个是，我发现如果不把所有的 $\Delta x^2+\Delta y^2$ 取成一个定值，画出来的图就很不好看。参考代码如下：

clc,clear,close all
x_0=-3:0.2:3;
y_0=-3:0.2:3;
[x,y]=meshgrid(x_0,y_0);
d=sqrt(1+(1-y.^2).^2);
u=ones(size(x));
v=(1-y.^2);
quiver(x,y,u,v);
xlim([-3,3])
ylim([-3,3])

所以我猜测，确实所有的 $\Delta x^2+\Delta y^2$ 都要取成一个定值。而且测试表明，这个定值不影响作图的效果。所以不妨设这个定值为1.

结合式子 $(1),(2)$，可以得到：

$$\Delta y=\frac{y^{\prime }}{\sqrt{1+y^{\prime 2}}},\Delta x=\frac{1}{\sqrt{1+y^{\prime 2}}}$$

这就是第5,6,7行代码的由来

---

2022年4月27日23:05:45