Matlab求解微分方程并画图

matlab求解微分方程并画图~

由于方程比较复杂，解析解不能用初等函数表示
只是要获得图像的话，用数值计算的办法可能更方便

fun=@(t,y) 1.44*(10^9)*(1-y).*exp(-109170./(8.314*t));
[T,Y] = ode23t(fun,500:600,0);
DY=fun(T,Y);
AX=plotyy(T,Y,T,DY);
set(get(AX(1),'Xlabel'),'String','T');
set(get(AX(1),'Ylabel'),'String','\alpha');
set(get(AX(2),'Ylabel'),'String','d\alpha/dT');

T是自变量，Y是变量也就是alpha
DY是，Y对T的导数
得到得到结果如下图



从图像看出，在T=500时，alpha=0
大概在T=600时，d(alpha)/dT趋向于0，alpha趋向于1不变
但是由于是数值解，在T不断增加的时候，d(alpha)/dT会在0附近振荡
所以T的取值不宜取得太大，这里取500到600之间
而采用ode23t函数，减少振荡

实际上，当T增大到一定值的时候，d(alpha)/dT趋向于0，
函数趋向于保持恒定值，所以后续的振荡是不合理的
取一定的区间如上图，已经可以很好地得到函数的变化趋势了

由于方程比较复杂，解析解不能用初等函数表示
只是要获得图像的话，用数值计算的办法可能更方便

fun=@(t,y) 1.44*(10^9)*(1-y).*exp(-109170./(8.314*t));
[T,Y] = ode23t(fun,500:600,0);
DY=fun(T,Y);
AX=plotyy(T,Y,T,DY);
set(get(AX(1),'Xlabel'),'String','T');
set(get(AX(1),'Ylabel'),'String','\alpha');
set(get(AX(2),'Ylabel'),'String','d\alpha/dT');

T是自变量，Y是变量也就是alpha
DY是，Y对T的导数
得到得到结果如下图



从图像看出，在T=500时，alpha=0
大概在T=600时，d(alpha)/dT趋向于0，alpha趋向于1不变
但是由于是数值解，在T不断增加的时候，d(alpha)/dT会在0附近振荡
所以T的取值不宜取得太大，这里取500到600之间
而采用ode23t函数，减少振荡

实际上，当T增大到一定值的时候，d(alpha)/dT趋向于0，
函数趋向于保持恒定值，所以后续的振荡是不合理的
取一定的区间如上图，已经可以很好地得到函数的变化趋势了

看标题以为你要求微分方程呐，结果是画dR/dr vs. r

% 画出图中的公式
% 定义微分方程函数
dRdr = @(r) 0.89 ./ r .* exp(-(log(r) + 0.84).^2 / 0.086);
% 在(0, 10]上画图
r = linspace(0.01, 3, 500);
dR = dRdr(r);
figure(10);
plot(r, dR, '^b', 'MarkerFaceColor', 'blue', 'DisplayName', 'dR/dr vs r');
xlabel('r', 'FontSize', 16);
ylabel('dR/dr', 'FontSize', 16);
legend('show');

作图的结果是

-----------------如果要求解微分方程的话-------------------------

这是一个一阶线性微分方程, 可以用龙格库塔求解. 首先用一个.m文件

定义图中的微分方程. 然后再用matlab的ode45函数求解这个方程.

举个例子, 建立一个solveFcn.m的文件如下

% 调用MATLAB的`ode45`函数实现求解. 具体说明在
% MATLAB命令行中输入`doc ode45`查询. 
% 简单来说, 需要两步. 1) 微分方程定义; 2) ode45参数设置和调用.
% 具体如下.
% 设置ode45参数, 并调用求解微分方程`DR = odefunc(r, R)`. 该方程定义
% 见后一个函数.
function Solution = solveFcn(initial_value)
% initial_value = [r0, R0];
r0 = initial_value(1);
R0 = initial_value(2);
% 由于提供的初值[r0 R0] = [0.43 0.5]; 或 [0.55 0.9]. 如果求解区间从(0, 10]
% 那么需要分成(0 r], [r, 10]两部分求解.
% [r0, 10]上求解
sol1 = ode45(@odefunc, [r0 10], R0);
% (0, r0]上求解
sol2 = ode45(@odefunc, [r0 0.01], R0); % log(r)在0处inf, 所以用0.01代替.
% sol结构体是微分方程的解, 包含`r` - 自变量, `R` - 要求解的函数R(r)
% 的数值. 其它r对应的R值, 可以用 `R = deval(sol, r);`计算. 其结果与
% sol中存储的r, R(r)相对精度一致.
Solution.sol1 = sol1;
Solution.sol2 = sol2;
end
% 定义微分方程. 用DR表示dR/dr. 一阶线性微分方程的一般形式为DR = f(r, R)
% 尽管题目微分方程右边未出现R, 此变量也不可省略.
function DR = odefunc(r, R)
DR = 0.89 ./ r .* exp(-(log(r) + 0.84).^2 / 0.086);
end

这个函数的作用就是求解最后几行定义的 dR/dr = f(r, R)这个微分方程. 解存在sol结构体中.(这个例子是在sol1和sol2里面)

怎么用这个sol解具体算每个r处的R呢, 需要用到deval. 可以再建立一个main_0.m文件, 内容为

% 调用函数sovleFcn. 求解微分方程
% 初值 r= 0.43, R = 0.5
S1 = solveFcn( [0.43 0.5] );
% 绘图
r1 = linspace(0.01, 0.43, 50); 
R1 = deval(S1.sol2, r1);
r2 = linspace(0.43, 2, 50);
R2 = deval(S1.sol1, r2);
r = [r1, r2];  % 在 0.01 - 2上的解
R = [R1, R2]; % 在 0.01 - 2上的解
figure(1); hold on;
plot(r, R, 'b-', 'linewidth',2, 'DisplayName', 'r=0.43, R=0.5');
% 初值 r=0.55, R = 0.9
S2 = solveFcn( [0.55 0.9] );
% 绘图
r1 = linspace(0.01, 0.55, 50); 
R1 = deval(S2.sol2, r1);
r2 = linspace(0.55, 2, 50); 
R2 = deval(S2.sol1, r2);
r = [r1, r2];  % 在 0.01 - 2上的解
R = [R1, R2]; % 在 0.01 - 2上的解
figure(1); hold on;
plot(r, R, 'r-', 'linewidth',2, 'DisplayName', 'r=0.55, R=0.9');
legend('show')
xlabel('r', 'FontSize', 16);
ylabel('\int_0^r dR', 'FontSize', 16)

那么这个微分方程在不同的初值条件下的解, 就可以画成下图

---

五原县13397665103： 如何用matlab求解微分方程并画图 - ？
夹很希捷： ^如何用matlab求解程并画图 dsolve()求解微分方程后,使用ezplot()作图,参考代码:1 2>> y = dsolve('Dy=y*cos(x)/(1+2*y^2)','y(0)=1','x'); >> ezplot(y,[-10,10]) % [-10,10]为x区间

五原县13397665103： 用matlab把微分方程解出来并画出图像!我要代码!大神!谢谢大神! - ？
夹很希捷： 此微分方程组可以用数值的方法求得其数值解.用ode45()函数,得到微分方程组解的图形. 执行代码和微分方程组定义函数func(t,y)

五原县13397665103： MATLAB求解微分方程并且画图,x为输入,F为输出 - ？
夹很希捷： 求解微分方程的函数是dsolve().如要画微分方程的解函数,可以用plot()函数

五原县13397665103： matlab解微分方程组并画图 - ？
夹很希捷： http://zhidao.baidu.com/question/205690281.html function dx=appollo(t,x) mu=1/82.45; mustar=1-mu; r1=sqrt((x(1)+mu)^2+x(3)^2); r2=sqrt((x(1)-mustar)^2+x(3)^2); dx=[x(2)2*x(4)+x(1)-mustar*(x(1)+mu)/r1^3-mu*(x(1)-mustar)/r2^3 x(4)-2*x(2)+x(3)-...

五原县13397665103： 用matlab(画图)求解微分方程(代码) sos啊!!! - ？
夹很希捷： matlab code:clear all clc f=@(t,y)([-y(2)+y(1)*(1-y(1)^2-y(2)^2); y(1)+y(2)*(1-y(1)^2-y(2)^2)]); [t y]=ode45(f,0:0.1:10,[5 5]); plot(t,y(:,1),t,y(:,2),y(:,2),y(:,1)) legend('x-t','y-t','x-y')

五原县13397665103： matlab解二阶微分方程并作图 - ？
夹很希捷： 假设a=b=c=1 dsolve('D2x+x-sin(t)')ans =(3*sin(t))/8 - sin(3*t)/8 - cos(t)*(t/2 - sin(2*t)/4) + C2*cos(t) + C3*sin(t)画图的有初值没有初值 C2 C3无法确定,无法画图

五原县13397665103： matlab 解微分方程并画出图 - ？
夹很希捷： 先定义函数,再用ode15s或ode45解.建一.m文件,在里面输入:function dy=xiao(x,y) dy=zeros(1,1); dy(1)=5.42*exp(-5.42)-2.5*y(1); 保存为xiao.m 再在command window界面输入:[x,y]=ode15s(@xiao,[0 10],[0]); plot(x,y) 即可.[0]对应于y的初值,[0 10]为一范围,都可自己设.ode15s改成ode45也可.

五原县13397665103： 怎么用MATLAB求解微分方程组并画出解函数图 - ？
夹很希捷： 方法如下:function dz=caonim(t,z) %保存为 caonim.m 文件 dz=[-z(1)^3-z(2);z(1)-z(2)^3];------------------------------------ 运行 [T,Z]=ode45('caonim',[0 30],[1,0.5]); plot(T,Z) legend('x-t','y-t')------------------------

五原县13397665103： 求matlab解一阶微分方程组,并画相图.dx/dt=y; dy/dt= - 2.6*y^3+0.6*y - 0.1*x - ？
夹很希捷：[答案] 你没给初值,我随便设一个[1 1] 新建一个脚本文件,输入: clear all clc f=@(t,x)([x(2);-2.6*x(2)^3+0.6*x(2)-0.1*x(1)]); [t,x]=ode45(f,[0 1],[1 1]); plot(t,x(:,1),t,x(:,2)) xlabel('t') legend('x','y')

五原县13397665103： 用MATLAB求微分方程,并画出图形解y=dsolve('D3y=D2y - Dy - y+t^2','y(0)=0','Dy(0)=1','D2y(0)= - 1')这样解了之后怎么画出图形啊! - ？
夹很希捷：[答案] y=dsolve('D3y=D2y-Dy-y+t^2','y(0)=0','Dy(0)=1','D2y(0)=-1'); tt=linspace(0,10,1001); %表示t从0到10 for i=1:1001 t=tt(i); yy(i)=eval(y); end plot(tt,yy)