动手练习

毫无章法系列（ode45用法真的别有洞天，解二阶微分方程组根本不在话下））
solve、dsolve
符号向量函数
ode函数
- 一阶微分方程
- 二阶微分方程组要转化为一阶微分方程组
- - 1
  - 2
  - 3
  - 4
  - 5
- 计算和扩展解结构体
vpa，simplify，subs
边界值问题的微分方程

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

solve、dsolve

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

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) 的表达式

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')

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)

符号向量函数

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

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)]

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等价，正向旋转

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

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

1

相应的函数句柄是微分方程组，初始值向量是多元组
总之要将求解的微分方程化成
y 1 ′ = f ( y 1 ) y_1'=f(y1) y1′=f(y1)
y 2 ′ = f ( y 2 ) y_2'=f(y2) y2′=f(y2)
y 1 ( 0 ) = 初始值 1 y_1(0)=初始值1 y1(0)=初始值1
y 2 ( 0 ) = 初始值 2 y_2(0)=初始值2 y2(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)

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

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

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('两者之间的距离')

计算和扩展解结构体

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
官方用法解释
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)]; endfunction res=dropbc(ya,yb);
res=[ya(1);yb(1)]; endfunction yinit=guess(x);
yinit=[sqrt(1-x.^2);-x./(0.1+sqrt(1-x.^2))]; end

微分方程解析解+数值解相关推荐

数学实验2：微分方程的数值解与解析解
文章目录一.实验目的和要求二.实验内容需要word文件请访问 http://daxs.top 站内搜索实验名称或者实验内容访问文章并且下载附件即可. 一.实验目的和要求熟练掌握掌握微分方程的数 ...
matlab中使用ode方法解范德波尔微分方程的数值解
微分方程的解析解要求比较严苛,只有在特定的条件下才能写出解析解表达式,而在现实的科研问题当中,绝大多数情况我们会采用数值解(numeric solution)的方法来求解微分方程.这个时候就要用到od ...
Matlab求微分方程的数值解
注:首先计算微分方程的解析解,如果发现没有解析解,再用数值解一.Matlab中求微分方程的数值解函数 [x,y]=solver('f',ts,x0,options) 1)x代表自变量 2)y代表函数 ...
matlab求解微分方程的数值解
简介:前面介绍了微分方程的解析解方法,同时也指出很多非线性微分方程是不存在解析解法的,需要使用数值解法对之进行研究.本文着重讨论基于 MATLAB/Simulink语言的各类微分方程的数值解方法. ...
计算物理中matlab处理微分方程解析解和欧拉法数值解的算法演示
我先来看一个问题的引入: 我们根据题目给出的微分方程编写matlab求解代码如下: syms cd m g v(t); eqn = diff(v,t) == g - cd/m*v^2; vt = ds ...
matlab--常微分方程的数值解(ODE-s)
用ODE(ordinary differential equations)23和ODE45求解微分方程 ODE(ordinary differential equations)23 [t,y]=ode ...
常微分方程数值解matlab欧拉,MATLAB题，用到欧拉公式求微分方程的数值解
%欧拉法解一阶常微分方程 %例子dy/h=-y+x+1 %f=inline('-y+x+1','x','y'); %微分方程的右边项 f = inline('x-2*y','x','y'); y0 ...
Matlab答疑五：使用微分定义求解微分方程的数值解
1.题目解微分方程 dydt=sin(y)+t,其中t=0时y=0,并绘图. 说明,一般对dydt的求解方法为:y(t+dt)=y(t)+dydt(t)*dt 2.方法除了题目给出方法:使用定义求 ...
matlab求解微分方程解析解
解析解 1求解微分方程方程组用函数"dsolve" dsolve('方程1','方程2'-'方程n','初始条件','自变量') clc,clear %% %求解du/dt=i+u ...

微分方程解析解+数值解