1 MATLAB常微分方程求解

1 找函数的解析解

1.1使用desolve函数，求解一阶微分方程

y1=dsolve('Dy==5'); 
y2=dsolve('Dy==x','x'); 
[y5,y6]=dsolve('Dx==y+x','Dy==2*x'); 
[y7,y8]=dsolve('Dx==y+x','Dy==2*x','x(0)==0','y(0)==1') ;
y9=dsolve('Dy==-2*y+2*x^2+2*x','y(0)==1','x') ;

%引号里的字符串都可以直接用变量代替
eqn10='Dy = y*x';
y10=dsolve(eqn10,'y(1)=1','x');

%画图方法
y10;
x = linspace(0,1,20);
z = eval(vectorize(y10));
plot(x,z)

结果

y1 = C1 + 5*t 
y2 = x^2/2 + C2
y5 = C8*exp(2*t) - (C7*exp(-t))/2
y6 = C7*exp(-t) + C8*exp(2*t)
y7 = exp(2*t)/3 - exp(-t)/3 
y8 = (2*exp(-t))/3 + exp(2*t)/3 
y9 = exp(-2*x) + x^2 
y10 = exp(-1/2)*exp(x^2/2) 

1.2使用desolve函数，求解二阶微分方程

eqn = 'D2y + 8*Dy + 2*y = cos(x)'; %D2y表示y''   Dy表示y'
inits = 'y(0) = 0,Dy(0) = 1';
y = dsolve(eqn,inits,'x'); 
x = linspace(0,10,300);
z = eval(vectorize(y)); 
plot(x,z); 

结果

y = (14^(1/2)*exp(4*x - 14^(1/2)*x)*exp(x*(14^(1/2) - 4))*(sin(x) - cos(x)*(14^(1/2) - 4)))/(28*((14^(1/2) - 4)^2 + 1)) - (14^(1/2)*exp(4*x + 14^(1/2)*x)*exp(-x*(14^(1/2) + 4))*(sin(x) + cos(x)*(14^(1/2) + 4)))/(28*((14^(1/2) + 4)^2 + 1)) - (14^(1/2)*exp(-x*(14^(1/2) + 4))*(7*14^(1/2) + 27))/(28*(8*14^(1/2) + 31)) - (14^(1/2)*exp(x*(14^(1/2) - 4))*(393*14^(1/2) + 1531))/(28*(8*14^(1/2) - 31)*(8*14^(1/2) + 31)^2)

图像

1.3 使用desolve函数，求解方程组

inits='x(0)=1,y(0)=2,z(0)=3';
[x,y,z]=dsolve('Dx=x+2*y-z','Dy=x+z','Dz=4*x-4*y+5*z',inits);
t = linspace(0,1,50);
xx = eval(vectorize(x))
yy = eval(vectorize(y))
zz = eval(vectorize(z))
plot(t,xx,'b',t,yy,'y',t,zz,'g') 

解得：

x = 6*exp(2*t) - (5*exp(3*t))/2 - (5*exp(t))/2  
y = (5*exp(3*t))/2 - 3*exp(2*t) + (5*exp(t))/2  
z = 10*exp(3*t) - 12*exp(2*t) + 5*exp(t)

2 找函数的数值解
2.1 ode45函数 百度百科传送门 https://baike.baidu.com/item/ode45/6674723?fr=aladdin

2.2 inline函数解决一阶微分方程

代码

f=inline('x*y^2+y');
[x,y]=ode45(f,[0 .5],1);
plot(x,y);

结果

增大步长（不设置步长的话，我也不知道初始步长是多少）
代码

xvalues=0:0.1:0.5
f=inline('x*y^2+y');
[x,y]=ode45(f,xvalues,1)
plot(x,y);

图像

2.3 误差的控制
我们在使用ode45算法的时候，运算结果是有误差的。RelTol表示绝对误差AbsTol表示相对误差。ek表示k点的误差，yk表示k点的函数值。存在以下关系

当yk很大的时候，我们的误差就会很大，就会得到错误的结果。
实验案例如下。
代码

xvalues=0:0.001:1
f=inline('x*y^2+y');
[x,y]=ode45(f,xvalues,1)
plot(x,y);

图像

可以看出， 如果不控制RelTol的话，我们计算的最大值是

>> max(y)

ans =

   1.0033e+03

为了得到精确的结果，我们控制一下相对误差，让结果更加准确
代码

xvalues=0:0.001:1
options=odeset('RelTol',1e-10);
f=inline('x*y^2+y');
[x,y]=ode45(f,xvalues,1,options)
plot(x,y);

图像

查看最大值，这个值就比较准确了

>> max(y)

ans =

   2.4695e+07

2.4 把函数代码放到.m文件里。
还是这个式子
firstode.m

function yprime = firstode(x,y);
% FIRSTODE: Computes yprime = x*yˆ2+y
yprime = x*y^2 + y;

在使用ode45时候，我们要加一个 @ 符号就好了
test.m

xspan = [0,.5];
y0 = 1;
[x,y]=ode23(@firstode,xspan,y0);
plot(x,y)

2.5 求解方程组


代码
lorenz.m

function xprime = lorenz(t,x);
%LORENZ: Computes the derivatives involved in solving the
%Lorenz equations.
sig=10;
beta=8/3;
rho=28;
xprime=[-sig*x(1) + sig*x(2); rho*x(1) - x(2) - x(1)*x(3); -beta*x(3) + x(1)*x(2)];

运行测试
test.m

x0=[-8 8 27];
tspan=[0,20];
[t,x]=ode45(@lorenz,tspan,x0)
plot(x(:,1),x(:,3)); %x(:,1)表示X  x(:,2)表示Y   x(:,3)表示Z

运行结果

test.m

 x0=[-8 8 27];
tspan=[0,20];
[t,x] = ode45(@lorenz,tspan,x0)
subplot(3,1,1);
plot(t,x(:,1));
subplot(3,1,2);
plot(t,x(:,2));
subplot(3,1,3);
plot(t,x(:,3));

运行结果

2.6 传递参数
lorenz1.m

function xprime = lorenz1(t,x,p);
%LORENZ: Computes the derivatives involved in solving the
%Lorenz equations.
sig=p(1); beta=p(2); rho=p(3);
xprime=[-sig*x(1) + sig*x(2); rho*x(1) - x(2) - x(1)*x(3); -beta*x(3) + x(1)*x(2)];

test.m

x0=[-8 8 27];
tspan=[0,20];
p=[10 8/3 28];
[t,x]=ode45(@lorenz1,tspan,x0,[],p);
plot(x(:,1),x(:,3))

应为传进去的参数和上次那个参数一模一样，所以这里就不贴图了。

2.7 二阶常微分方程求解

对于阶常微分方程的求解，可以通过增加变量个数的方式：Z=Y'，然后增加方程个数，写成方程组的形式。

实现代码
F.m

function xprime = F(t,x);
%LORENZ: Computes the derivatives involved in solving the
%Lorenz equations. 
xprime=[ 1 ; x(3) ; -8*x(3)-2*x(2)+cos( x(1) ) ]; 

test.m

x0=[0 0 1];
tspan=[0,10];
[t,x] = ode45(@F,tspan,x0)
plot(x(:,1),x(:,2)); %x(:,1)表示X  x(:,2)表示Y   x(:,3)表示Z

图像

3 拉普拉斯转换
3.1拉普拉斯变换
示例代码

>>syms t;
>>laplace(tˆ2)

3.2拉普拉斯逆变换

>>syms s;
>>ilaplace(1/(1+s))

4 边界值问题(我们不知道所有的初始值，但是知道边界值)
对于以下方程组：

代码
bc.m

function res=bc(L,R)
%BC: Evaluates the residue of the boundary condition
res=[L(1)-3;R(1)-10];  %L = 1    R= 10要写成这种形式。

bvpexample.m

function yprime = bvpexample(t,y)
%BVPEXAMPLE: Differential equation for boundary value
%problem example.
yprime=[y(2); -2*y(1)+3*y(2)];

test.m

sol=bvpinit(linspace(0,1,25),[0 1]);
sol=bvp4c(@bvpexample,@bc,sol); 
plot( sol.x,sol.y(1,:) )

没学的部分1（资料不全，没Example3.1没办法先不学了）

没学的部分2 这一部分讲泰勒展开，高数老师没讲过，主要是一些数学原理，时间原因，先不学习。