对于微分方程
y ′ = 1 − y 2 y'=1-y^2 y′=1−y2
画出他的相图,其中 x ∈ [ − 3 , 3 ] x\in[-3,3] x∈[−3,3] , y ∈ [ − 3.3 ] y\in [-3.3] y∈[−3.3],步长为 h = 0.2 h=0.2 h=0.2,箭头长为 c = 0.01 c=0.01 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 ] × [ − 3 , 3 ] [-3,3]\times[-3,3] [−3,3]×[−3,3] 上等距取点,也就是第二三四行代码。
那么箭头的方向和大小如何计算呢?首先要知道箭头的方向是什么含义。箭头的方向表示曲线在该点处的切线方向,而一点处的切线方向,可以用 Δ x , Δ y \Delta x,\Delta y Δx,Δy 来表示。
而我们要算的 u v,就是
Δ
x
,
Δ
y
\Delta x,\Delta y
Δx,Δy。
由微积分的知识我们可以知道。
Δ y Δ x = y ′ (1) \frac{\Delta y}{\Delta x}=y'\tag{1} ΔxΔy=y′(1)
可以想象到在该切线方向上,有任意多满足 ( 1 ) (1) (1) 的 ( Δ x , Δ y ) (\Delta x, \Delta y) (Δx,Δy) ,但我们这里考虑单位向量,即:
Δ x 2 + Δ y 2 = 1 (2) \Delta x^2+\Delta y^2=1\tag{2} Δx2+Δy2=1(2)
至于我为什么要考虑单位向量呢?原因有二,一个是总要有一个标准来确定 Δ x , Δ y \Delta x,\Delta y Δx,Δy,选啥向量都一样,还不如选个单位向量感觉比较酷。另一个是,我发现如果不把所有的 Δ x 2 + Δ y 2 \Delta x^2+\Delta y^2 Δx2+Δy2 取成一个定值,画出来的图就很不好看。参考代码如下:
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])
所以我猜测,确实所有的 Δ x 2 + Δ y 2 \Delta x^2+\Delta y^2 Δx2+Δy2 都要取成一个定值。而且测试表明,这个定值不影响作图的效果。所以不妨设这个定值为1.
结合式子 ( 1 ) , ( 2 ) (1),(2) (1),(2),可以得到:
Δ y = y ′ 1 + y ′ 2 , Δ x = 1 1 + y ′ 2 \Delta y= \frac{y'}{\sqrt{1+y'^2}},\Delta x=\frac{1}{\sqrt{1+y'^2}} Δy=1+y′2y′,Δx=1+y′21
这就是第5,6,7行代码的由来
2022年4月27日23:05:45