如何用matlab求解二阶微分方程，以及程序实例

如何用matlab求解二阶微分喊冲态方程，对于一般的微分方程（组）可以用dsolve（）函数求得其解析解或数判绝值解，对于较复杂的微分方程（组）可以用ode45（）函郑源数求得其数值解。

例如：微分方程，用dsolve和ode45计算t=0.1,0.2时y的值。

y = dsolve('D2y+0.5*Dy+2*y =0','Dy(0)=0,y(0)=1')

t=0.1y1=eval(y)

t=0.2y2=eval(y)

==============================================================

t0=[10]tspan=[0,20]

[t1,y1]=ode45(@ ode_fun,tspan,t0)；

disp('t      y      dy')

for i=1:3

   t=double(t1(i))y=double(y1(i,1))dy=double(y1(i,2))

   disp([num2str(t),'      ',num2str(y),'      ',num2str(dy)]) 

end

function f = ode_fun(t,y)   %自定义的微分方程函数，文件名，ode_fun.m

f=[y(2)-0.5*y(2)-2*y(1)]

end

==============================================================

运行结果

主程序：

clear

%1 已知参数

alpha_0=0phi_0=0

theta_1=[4590135]*pi/180

theta_3=[5282112]*pi/180

d=50

%2 计算各杆长度

[l,m,n]=link_design(theta_1,theta_3,alpha_0,phi_0)

b=l*d/n

c=d*m/n

a=d/n

%3 输出计算结果

disp('计算结果1：各干相对长度')

disp('            ')

fprintf('连杆相对长度 l=%3.2f  \n',l)

fprintf('摇杆相对长度 m=%3.2f  \n',m)

fprintf('机架相对长度 察基n=%3.2f  \n'凯销,n)

disp('            ')

disp('计算结果2：各干长度')

disp('            ')

fprintf('曲柄长度 a=%3.2f \n',a)

fprintf('连杆长度 b=%3.2f \n',b)

fprintf('摇杆长度 c=%3.2f \n',c)

fprintf('机架长度 d=%3.2f \n',d)

disp('            ')

子程序（link_design,m）

function   [l,m,n]=link_design(theta_1,theta_3,alpha_0,phi_0)

%计算线性方程组系数矩阵A

A=[cos(theta_3(1)+phi_0),cos(theta_3(1)+phi_0-theta_1(1)-alpha_0),1

    cos(theta_3(2)+phi_0),cos(theta_3(2)+phi_0-theta_1(2)-alpha_0),1

    cos(theta_3(3)+phi_0),cos(theta_3(3)+phi_0-theta_1(3)-alpha_0),1]

% 计算线性方程组系数矩阵B

B=[cos(theta_1(1)+alpha_0)cos(theta_1(2)+alpha_0)cos(theta_1(3)+alpha_0)]

% 求解线性方程组

p=A\B

% 计算相对杆长l,m,n

p0=p(1)p1=p(2)p2=p(3)

m=p0

n=-m/p1

l=sqrt(m*m+n*n+1-p2*2*n)

结果

计算结果1：各干相对长度

            

连杆相对长度 l=2.07  

摇杆相对长度 m=1.49  

机架相对长度 n=1.81  

            

计算结果2：各败孙谨干长度

            

曲柄长度 a=27.63 

连杆长度 b=57.24 

摇杆长度 c=41.11 

机架长度 d=50.00

clear

clc

a = 0.95

k = [510134311131081674]

k = -k % 模拟退火算法是求解最小值，故取负数

d = [2518325104117146]

restriction = 46

num = 12

sol_new = ones(1,num)% 生成初始解

E_current = infE_best = inf

% E_current是当前解对应的目标函数值（即背包中物品总价值销猜穗）；

% E_new是新解的目标函兆告数值；

% E_best是最优解的

sol_current = sol_newsol_best = sol_new

t0=97tf=3t=t0

p=1

while t>=tf

for r=1:100

%产生随机扰动

tmp=ceil(rand.*num)

sol_new(1,tmp)=~sol_new(1,tmp)

%检查是否满足约束

while 1

q=(sol_new*d <= restriction)

if ~q

p=~p %实现交错着逆转头尾的第一个1

tmp=find(sol_new==1)

if p

sol_new(1,tmp)=0

else

sol_new(1,tmp(end))=0

end

else

break

end

end

% 计算背包中的物品亏卜价值

E_new=sol_new*k

if E_new<E_current

E_current=E_new

sol_current=sol_new

if E_new<E_best

% 把冷却过程中最好的解保存下来

E_best=E_new

sol_best=sol_new

end

else

if rand<exp(-(E_new-E_current)./t)

E_current=E_new

sol_current=sol_new

else

sol_new=sol_current

end

end

end

t=t.*a

end

disp('最优解为：')

sol_best

disp('物品总价值等于：')

val=-E_best

disp(val)

disp('背包中物品重量是：')

disp(sol_best * d)

---