matlab解微分方程

1.匿名函数

1.1创建

f=@(变量) 表达式;

f=@(x1,x2) x1.^2+x2;

1.2 求解

x1为2 3 4 5；x2为3 4 5 6的情况下求解函数f的值

f=@(x1,x2) x1.^2+x2;
y=f(2:5,3:6);
subplot(121);%选择子图位置
plot(y)%画图

2.一阶微分方程

用“dsolve”

2.1例 y.-y=0

step1: 申明自变量和因变量

• syms y(x)

step2：编程

clear all;
syms y(x)
ode=diff(y,x)-y==0;				%注意是==
cond=[];
dsolve(ode,cond);

得到：

ans =
 
C1*exp(x)

2.2例 y.-y=x

step1: 申明自变量和因变量

• syms y(x)

step2：编程

clear all;
syms y(x)
ode=diff(y,x)-y==x;				%注意是==
cond=[];
dsolve(ode,cond);

得到

ans =
 
C1*exp(x) - x - 1

2.3例dy/dx=sinx

step1.写函数文件

function dy = ST1(x,y)
dy= sin(x);
end

step2.主函数

clear all;
[X,Y]=ode45('ST1',[2*pi],[1]);
plot(X,Y(:,1))

3.二阶微分方程

相当于定义了一个新向量y，然后列匿名函数，方程的左边都是一阶导数，如y(1).;y(2).等等，右边直接打公式。

3.1例，x… = -x，x初值0，x.初值2

step1：改写

• y(1)=x；
• y(2)=x.

step2：得y的导数

• y(1).=x.=y(2)
• y(2).=x…=-x=-y(1)

step3：编程

clear;
F=@(t,y) [0,1;-1,0]*y;
ode45(F,[0,10],[0,2])%[0,10]表示t的范围，0，2表示积分初值

3.2例 ax…+bsin(x)=0

step1：改写

• y(1)=x;
• y(2)=x.;

step2：得y的导数

• y(1).=x.=y(2) ;
• y(2).=x…=-bsin(x)/a;

step3：编程

clear;
F=@(t,y)[y(2);-9.8*sin(y(1))];
[t,y]=ode45(F,[0,20],[9.5*pi/10,0]);%注意这里加了[t,y]
plot(y(:,1),y(:,2));%这样才能输出y的值

3.3例洛伦兹方程

step1：改写

• y(1)=x
• y(2)=y
• y(3)=z

step2：得y的导数

• y(1).=-ay(1)+ay(2);
• y(2).=b*y(1)-y(2)-y(1)*y(3);
• y(3).=y(1)y(2)-cy(3);

step3：编程

clear;
a=10;
b=35;
c=8/3;
F=@(t,y)[-a*y(1)+a*y(2);
    b*y(1)-y(2)-y(1)*y(3);
    y(1)*y(2)-c*y(3)];
[t,y]=ode45(F,[0,20],[-0.01,0.0,0.001]);
plot3(y(:,1),y(:,2),y(:,3));

3.4 耦合

step1：改写

• y(1)=r
• y(2)=r.
• y(3)=xita
• y(4)=xita.

step2：得y的导数

• y(1).=y(2)
• y(2).=-a/y(1)²+y(1) y(4)²
• y(3).=y(4)
• y(4).=-b*y(2)*y(4)/y(1)

step3：编程

clear;
a=5.965*6.67295;
b=2;
F=@(t,y)[y(2);
    -a/(y(1)^2)+y(1)*y(4)^2;
    y(4);
    -b*y(2)*y(4)/y(1)];
[t,y]=ode45(F,[0,1500],[63.71,0.0,0.001,0.79/63]);
[x,y]=pol2cart(y(:,3),y(:,1));%极坐标转化为直角坐标
plot(x,y)

4.混合求导的耦合微分方程

这个叫法肯定不太专业，主要是求导变量换了。

4.1例1

step1.写函数文件

function dy = ST1(x,y)
dy=zeros(2,1);%一定要写！！！！，声明列向量的维度！！！！
dy(1)= y(1)+y(2)/x;
dy(2)= sin(y(1))+x;
end

step2：主程序ode求解

clear all;
[X,Y]=ode45('ST1',[2*pi,3*pi],[1,1]);
plot(X,Y(:,1),X,Y(:,2))

4.2例

step1：改写

• y(1)=x
• y(2)=y

step2：得y的导数

• y(1).=r1 y(1) (1-y(1)/n1-s1*y(2)/n2)
• y(2).=r2 y(2) (1-y(2)/n2-s2*y(1)/n1)

step3：编程

clear;
r1=1;
r2=0.2;
n1=100;
n2=100;
s1=0.5;
s2=2;
F=@(t,y)[r1*y(1)*(1-(y(1)/n1)-s1*(y(2)/n2));
    r2*y(2)*(1-(y(2)/n2)-s2*(y(1)/n1))];
[t,y]=ode45(F,[0,20],[10,10]);
plot(t,y(:,1),t,y(:,2))