微分方程解析解+数值解

毫无章法系列（ode45用法真的别有洞天，解二阶微分方程组根本不在话下））

solve、dsolve

dsolve('方程1','方程2',...,'方程n','初始条件','自变量')
1

clc,clear;
syms m V rho g k;
%微分方程求解(解就是满足微分方程的 变量之间的关系式）
dsolve('Dy==5')  %dy/dt=5
dsolve('Dy==x','x') %dy/dx=x 没指定变量则默认为t
dsolve('D2y==1+Dy','y(0)==1','Dy(0)==0')%d2y/dt2=1+dy/dt,初始条件y(0)=1,dy(0)/dt=0
[y1,y2]=dsolve('Dx==y+x','Dy==2*x','x(0)==0','y(0)==1')%dx/dt=y+x,dy/dt=2*x,x(0)=0,y(0)=1,x=y5,y=y6

%提前定义了符号变量，就可以直接在函数括号里写表达式，不用加引号；
%没定义那么就放在引号里
%solve对一般方程的求解
syms x y a;
eq=x^2+2*x+1;
s1=solve(eq,x)
eq=a*x+2;
s2=solve(eq,a)
eq1=x+2*y-8; %解二元一次方程组，解是个结构体
eq2=3*x+5*y-4;
s3=solve(eq1,eq2,x,y);%只用s3承载结果则s3为向量组，可用[x,y]承接
s3=[s3.x,s3.y]

s1=solve(sin(x)==1/2)
s2=solve(x^3-1==0)

%ode用于求微分方程的数值解，实例见CSDN，重点在于
%微分方程的标准形式，总能找出多个变量 如F(y,y',y'',y''',…,t)=0
%y,y',y'',y'''就是变量，要设一个向量，做变量替换
%如向量x：  x(1)=y，x(2)=y'
%ode45 函数主要部分 要分别列出每个变量 对t求导的 d x(2) 的表达式


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

clc, clear, 
syms y(x)
y=dsolve(diff(y)==-2*y+2*x^2+2*x, y(0)==1)
disp('----------等同于-------------')
y=dsolve('Dy==-2*y+2*x^2+2*x','y(0)==1','x')

clc,clear,syms y(x) %定义符号函数y，自变量为x
dy=diff(y);
y=dsolve( diff(y,2)-2*diff(y)+y ==exp(x) , y(0)==1 ,dy(0)==-1 )
disp('----------等同于-------------')
y=dsolve('D2y-2*Dy+y==exp(x)','y(0)==1','Dy(0)==-1','x')


1
2
3
4
5
6
7
8
9
10
11
12
13

clc,clear,syms y(t)
u=exp(-t)*cos(t)
dy=diff(y);d2y=diff(y,2);d3y=diff(y,3);d4y=diff(y,4);
y=dsolve( d4y+10*d3y+35*d2y+50*dy+24*y==diff(u,2),...
    y(0)==0,dy(0)==-1,d2y(0)==1,d3y(0)==1 )
% y=dsolve(diff(y,4)+10*diff(y,3)+35*diff(y,2)+50*diff(y)+24*y==diff(u,2),...
%     y(0)==0,dy(0)==-1,d2y(0)==1,d3y(0)==1)


1
2
3
4
5
6
7
8
9

符号向量函数

建立符号向量函数参考博客

clc,clear;
%syms x(t) [3,1] %符号向量函数  不可行
syms x(t) y(t) z(t)
X=[x;y;z];
A=[3 -1 1;2 0 -1;1 -1 2];
[s1,s2,s3]=dsolve(diff(X)==A*X,X(0)==ones(3,1))
%[s1,s2,s3]=dsolve(diff(X)==A*X,X(0)==[1;1;1])
s=[simplify(s1);simplify(s2);simplify(s2)]

1
2
3
4
5
6
7
8
9

ode函数

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

• odefun 是函数句柄，可以是函数文件名，匿名函数句柄或内联函数名
• tspan 是区间 [t0 tf] 或者一系列散点[t0,t1,…,tf]
• y0 是初始值向量
• T 返回列向量的时间点
• Y返回对应T的求解列向量

一阶微分方程 ode45函数主体是 因变量
@(匿名函数的输入参数)

一阶微分方程

比较符号解和数值解，符号解的图hold on，数值解换一种图例接着画，于是数值解落在符号解上

clc, clear, close all, syms y(x)
y=dsolve(diff(y)==-2*y+2*x^2+2*x, y(0)==1)
dy=@(x,y)-2*y+2*x^2+2*x;
[sx, sy]=ode45(dy, [0,0.5], 1)
fplot(y,[0,0.5]), hold on    %符号解
plot(sx, sy, '*'); 
legend({'符号解','数值解'})   %图例
xlabel('$x$','Interpreter','Latex')
ylabel('$y$','Interpreter','Latex','Rotation',0)
%rotation代表旋转的意思，这里0和360等价，正向旋转
1
2
3
4
5
6
7
8
9
10

rotation代表旋转的意思，这里0和360等价，正向旋转

二阶微分方程组要转化为一阶微分方程组

1

相应的函数句柄是微分方程组，初始值向量是多元组
总之要将求解的微分方程化成
$y_1^{\prime }=f(y1)$
$y_2^{\prime }=f(y2)$
$y_1(0)=初始值1$
$y_2(0)=初始值2$

2

clc, clear, close all
dy=@(x,y)[y(2); sqrt(1+y(2)^2)/5/(1-x)];
[x,y]=ode45(dy,[0,0.999999],[0;0])
plot(x, y(:,1));
xlabel('$x$','Interpreter','Latex')
ylabel('$y$','Interpreter','Latex','Rotation',0)

1
2
3
4
5
6
7

3

function Testode45
tspan=[3.9 4.0]; %求解区间
y0=[8 2]; %初值
[t,x]=ode45(@odefun,tspan,y0);
plot(t,x(:,1),'-o',t,x(:,2),'-*')
legend('y1','y2')
title('y'' ''=-t*y + e^t*y'' +3sin2t')
xlabel('t')
ylabel('y')
function y=odefun(t,x)
y=zeros(2,1); % 列向量
y(1)=x(2);
y(2)=-t*x(1)+exp(t)*x(2)+3*sin(2*t); %常微分方程公式
end
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

4

clear,clc
odefun=@(t,y)[y(2); (1-y(1)^2)*y(2)-y(1)];
[t,y]=ode45(odefun,[0 20],[2;0]);
plot(t,y(:,1),'o',t,y(:,2),'o')
legend('y','y的一阶导')
% function dy=odefun(t,y) 还是不要写函数体的好
% dy=[y(2);
% (1-y(1)^2)*y(2)-y(1)];
% end
1
2
3
4
5
6
7
8
9

5

df=@(t,f,w)[……;……]
s1=ode45(@(t,x)df(t,x,w),[t0,tf],[0;0])
其实就是为了将w也当作一个变量（w是不同的定值）
看过下图就懂

 clear,clc,close all, w=20;
df=@(t,f)[w/sqrt((10+20*cos(t)-f(1))^2+...
    (20+15*sin(t)-f(2))^2)*(10+20*cos(t)-f(1))
    w/sqrt((10+20*cos(t)-f(1))^2+...
    (20+15*sin(t)-f(2))^2)*(20+15*sin(t)-f(2))];
t0=0; tf=5; %时间的终值是适当估计的
s1=ode45(df,[t0,tf],[0;0]) %求微分方程的数值解
d1=@(t)sqrt((deval(s1,t,1)-10-20*cos(t)).^2+...
    (deval(s1,t,2)-20-15*sin(t)).^2); %t时刻人和狗的距离
[tmin,fmin]=fminbnd(d1,0,tf) %求两者距离的最小值及对应的时间
t=0:0.1:tf; subplot(121), plot(t,d1(t))  %画出两者之间的距离
xlabel('$t$','Interpreter','latex'), ylabel('两者之间的距离')

% w=5; tf=60; 
% [t,s]=ode45(@(t,x)df(t,x,w),[t0,tf],[0;0]); %求微分方程的数值解
% d2=sqrt((s(:,1)-10-20*cos(t)).^2+(s(:,2)-20-15*sin(t)).^2); %计算两者距离
% subplot(122), plot(t,d2)  %画出两者之间的距离
% xlabel('$t$','Interpreter','latex'), ylabel('两者之间的距离')

disp('----因为要求w分别再不同取值下的结构，将w也当作变量-----------')

clc, clear, close all, w=20;
df=@(t,f,w)[w/sqrt((10+20*cos(t)-f(1))^2+...
    (20+15*sin(t)-f(2))^2)*(10+20*cos(t)-f(1))
    w/sqrt((10+20*cos(t)-f(1))^2+...
    (20+15*sin(t)-f(2))^2)*(20+15*sin(t)-f(2))];
t0=0; tf=5; %时间的终值是适当估计的
s1=ode45(@(t,x)df(t,x,w),[t0,tf],[0;0]) %求微分方程的数值解
d1=@(t)sqrt((deval(s1,t,1)-10-20*cos(t)).^2+...
    (deval(s1,t,2)-20-15*sin(t)).^2); %t时刻人和狗的距离
[tmin,fmin]=fminbnd(d1,0,tf) %求两者距离的最小值及对应的时间
t=0:0.1:tf; subplot(121), plot(t,d1(t))  %画出两者之间的距离
xlabel('$t$','Interpreter','latex'), ylabel('两者之间的距离')

w=5; tf=60; 
[t,s]=ode45(@(t,x)df(t,x,w),[t0,tf],[0;0]); %求微分方程的数值解
d2=sqrt((s(:,1)-10-20*cos(t)).^2+(s(:,2)-20-15*sin(t)).^2); %计算两者距离
subplot(122), plot(t,d2)  %画出两者之间的距离
xlabel('$t$','Interpreter','latex'), ylabel('两者之间的距离')
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

计算和扩展解结构体

sol = ode45(___) 返回一个结构体，您可以将该结构体与 deval 结合使用来计算区间 [t0 tf] 中任意点位置的解。您可以使用上述语法中的任何输入参数组合。
官方文档

通过ode45的解sol，可以通过deval得到任意x对应的数值解

扩展解

vpa，simplify，subs

syms x y z 
f=cos(x)^2-sin(x)^2
s1 = simplify(f) 
s1 = cos(2*x) 
% 文件名不要与matlab固有函数名重名！！！
clc,clear;
syms m V rho g k;

s=dsolve('m*D2s=m*g-rho*g*V-k*Ds','s(0)=0','Ds(0)=0')
%该微分方程只有一个变量s，或者说 微分方程的解就是s与t的关系式

s=subs(s,{m,V,rho,g,k},{239.46,0.2058,1035.71,9.8,0.6})
%  R = subs(S, old, new) 利用new的值代替符号！表达式！中old的值

s=vpa(s,10) 
%使用符号计算时得到的精确解会出现分数，可以用vpa转换为小数显示
%vpa作用对象可以是数值或者！表达式！（表达式中的数值精度）（有效数字）

%s=dsolve('m*DV==m*g-rho*g*V-K*V')

syms f(x)
a=1/99
x=sym(1/2)
y=vpa(x)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

边界值问题的微分方程

参考博文1
参考博文2
参考博文3
官方用法解释
exp6.15

function main_Optimalcpntrol
clc, clear, close all
% drop=@(x,y)[y(2);(y(1)-1)*(1+y(2)^2)^(3/2)];
% dropbc=@(ya,yb)[ya(1);yb(1)];
% guess=@(x)[sqrt(1-x.^2);-x./(0.1+sqrt(1-x.^2))];
% guess=@(x)[x.^2;2*x];
solinit=bvpinit(linspace(-1,1,20), @guess);
sol=bvp4c(@drop, @dropbc, solinit)
% fill(sol.x, sol.y(1,:), [0.7,0.7,0.7]), axis([-1,1,0,1])

xint = linspace(-1,1,20);
yint = deval(sol, xint);
plot(xint, yint(1,:));% y(1)是y，y(2）是y`
axis([-1,1,0,1])

xlabel('$x$','Interpreter', 'latex', 'FontSize',12)
ylabel('$h$','Interpreter', 'latex', 'Rotation', 0, 'FontSize',12)
end
function yprime=drop(x,y);
yprime=[y(2);(y(1)-1)*(1+y(2)^2)^(3/2)]; end
 
function res=dropbc(ya,yb);
res=[ya(1);yb(1)]; end
 
function yinit=guess(x);
yinit=[sqrt(1-x.^2);-x./(0.1+sqrt(1-x.^2))]; end

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27