EDA365电子论坛网

标题: 请问怎样用MATLAB求解微分方程组并画出解函数图？ [打印本页]

作者: greensmile 时间: 2020-4-20 10:27
标题: 请问怎样用MATLAB求解微分方程组并画出解函数图？
求解微分方程组，画出解函数图。

x'= -x^3-y, x(0)=1
y'= x-y^3, y(0)=0.5 0<t<30
请教高手指点给出程序与结果，谢谢

作者: yin123 时间: 2020-4-20 13:34
本帖最后由 yin123 于 2020-4-20 13:36 编辑

碰巧正好学习了微分方程解法
matlab内置的函数（数值解法）有ode系列，help一下就可以了
如果你上的是有关数字分析的课程，我猜老师的意思应该让你自己写代码
我分别尝试了Euler法梯形法四阶RK法代码如下：
以y‘=|sin2x、，y（0）=1为例：

编写一阶导数m文件：

%一阶常微分方程的m文件
function y1=myfun(x,y)
y1=abs(sin(2*x));

Euler法：

%向前差分Euler法解常微分方程
%fun：f（x）的一阶导数m文件
%x0,xt：自变量的初值和终值
%y0：f（x）在x0处的值
%h：步长
function [Eulerx,Eulery]=MyEuler(myfun,x0,xt,y0,h)
x=(x0:h:xt)';%定义自变量数组
y(1,: )=y0;%初始化因变量数组
for k = 1:length(x)-1
f=feval(myfun,x(k),y(k,: ));%计算每个迭代点的一阶导数值
f=f(: )';
y(k+1,: )=y(k,: )+h*f; %向前迭代
end
Eulerx=x;
Eulery=y;

梯形法：

%梯形法解常微分方程
%fun：f（x）的一阶导数m文件
%x0,xt：自变量的初值和终值
%y0：f（x）在x0处的值
%h：步长
function [Tixingx,Tixingy]=MyTixing(myfun,x0,xt,y0,h)
x=(x0:h:xt)';%定义自变量数组
y(1,: )=y0;%初始化因变量数组
Y(1,: )=y0;%中间量
for k=1:length(x)-1
f=feval(myfun,x(k),y(k,: ) );%计算每个迭代点的一阶导数值
f=f(: )';
Y(k+1,: )=y(k,: )+h*f; %y0(n+1)
F=feval(myfun,x(k+1),Y(k+1,: ));%计算每个迭代点的一阶导数值
F=F(: )';
y(k+1,: )=y(k,: )+0.5*h*(f+F); %y(n+1)
end
Tixingx=x;
Tixingy=y;

四阶RK法：

%四阶Runge-Kutta法解常微分方程
%fun：f（x）的一阶导数m文件
%x0,xt：自变量的初值和终值
%y0：f（x）在x0处的值
%h：步长
function [RK4x,RK4y]=RK4(myfun,x0,xt,y0,h)
x=(x0:h:xt)';%定义自变量数组[code]clear all;clc;
x0=0;xt=2*pi;%x0,xt：自变量的初值和终值
y0=1;%y0：f（x）在x0处的值
h=0.2;%h：步长
[Eulerx,Eulery]=MyEuler(@myfun,x0,xt,y0,h);%欧拉法
a=fliplr(Eulery');%求算f（22*pi)
y_Euler=a(1);
[Tixingx,Tixingy]=MyTixing(@myfun,x0,xt,y0,h);%梯形法
a=fliplr(Tixingy');%求算f（22*pi)
y_Tixing=a(1);
[RK4x,RK4y]=RK4(@myfun,x0,xt,y0,h);%四阶Runge-Kutta法法
a=fliplr(RK4y');%求算f（22*pi)
y_RK4=a(1);
[ode45x,ode45y]=ode45(@myfun,[x0:h:xt],[1]);%matlab内置函数
a=fliplr(ode45y');%求算f（22*pi)
y_ode45=a(1);
figure;
figure_FontSize=18; %修改字体
plot(Eulerx,Eulery,'ro',Tixingx,Tixingy,'g*',RK4x,RK4y,'bd',ode45x,ode45y,'ks','LineWidth',2);%做图
legend('Euler','Tixing','RK4','ode45');%图例
Y=[y_Euler,y_Tixing,y_RK4,y_ode45];
axis([x0 xt y0 max(Y)]); %定义坐标轴范围
label1={'方法';'Euler';'Tixing';'RK4';'ode45'};
label2={'y(2π）';y_Euler;y_Tixing;y_RK4;y_ode45};
% label=[label1,label2];
text(x0+0.1*(xt-x0),0.8*max(Y),label1); %显示y(2π）
text(x0+0.3*(xt-x0),0.8*max(Y),label2);

y(1,: )=y0;%初始化因变量数组
for k=1:length(x)-1
K1=feval(myfun,x(k),y(k,: ));%计算系数
K2=feval(myfun,x(k)+0.5*h,y(k,: )+0.5*K1');
K3=feval(myfun,x(k)+0.5*h,y(k,: )+0.5*K2');
K4=feval(myfun,x(k)+h,y(k,: )+h*K3');
y(k+1,: )=y(k,: )+(1/6)*h*(K1'+2*K2'+2*K3'+K4');%迭代
end
RK4x=x;
RK4y=y;

作者: CCxiaom 时间: 2020-4-20 13:39

楼主的图用Forcal绘制，代码：

!using["XSLSF"]; //使用命名空间XSLSF
//数组xArray存放x的值；ti为当前有效值的个数；tmax为ti对应的时间；tmin为起始时间。
xt(t:k:xArray,ti,tmax,tmin)=
{
k=(t-tmin)/(tmax-tmin)*ti-1,
if[k<0, k=0], if[k>=ti, k=ti-1],
xArray.getrai[k]
};
//数组yArray存放y的值；ti为当前有效值的个数；tmax为ti对应的时间；tmin为起始时间。
yt(t:k:yArray,ti,tmax,tmin)=
{
k=(t-tmin)/(tmax-tmin)*ti-1,
if[k<0, k=0], if[k>=ti, k=ti-1],
yArray.getrai[k]
};
//函数定义，用于计算微分方程组中各方程右端函数值，连分式法对微分方程组积分一步函数pbs1将调用该函数f。
//t为自变量，x,y为函数值，dx,dy为右端函数值（即微分方程的值）。
f(t,x,y,dx,dy)=
{
dx=-(x^3)-y,
dy=x-y^3
};
//用连分式法对微分方程组进行积分，获得理论值。
//t1,t2为积分的起点和终点。
//h,s为自动变量。
//模块变量：hf为函数f的句柄，要预先获得该句柄；Array为工作数组；step为积分步长；eps为积分精度。
积分(t1,t2:h,s:hf,Array,step,eps)=
{
s=ceil[(t2-t1)/step], //计算积分步数
h=(t2-t1)/s, //重新计算积分步长
{ pbs1[hf,t1,Array,h,eps], //对微分方程组积分一步
t1=t1+h //积分步长增加
}.until[abs(t1-t2)<h/2] //连续积分至t2
};
数据(:i,h,t1,t2,x,y,static,free:hf,Array,step,eps,max,xArray,yArray,ti,tmax,tmin)= //微分方程组积分获得数据
{
if[free,delete(Array),delete(xArray),delete(yArray),return(0)],
hf=HFor("f"), //模块变量hf保存函数f的句柄，预先用函数HFor获得该句柄
step=0.01,eps=1e-6, //积分步长step要合适，积分精度eps越小越精确，用于对微分方程组积分一步函数pbs1
Array=new[rtoi(real_s),rtoi(2*15)], //申请工作数组
max=1001,
xArray=new[rtoi(real_s),rtoi(max)], //申请工作数组
yArray=new[rtoi(real_s),rtoi(max)], //申请工作数组
Array.setra(0,1,0.5), //设置积分初值，通过模块变量Array传递，Array是一个数组
xArray.setra(0,1), //设置xArray的第一个值
yArray.setra(0,0.5), //设置yArray的第一个值
ti=1, h=step*3, tmin=0, tmax=0, t1=0, t2=h,
i=1,(i<max).while{
积分(&t1,t2), //从t1积分到t2
Array.getra(0,&x,&y), //从数组Array获得t2时的积分值
xArray.setra(i,x), yArray.setra(i,y), //将积分值保存到数组
ti=i+1, tmax=t1, t2=t2+h,
i++
}
};
//绘制函数图形
(::hxt,hyt)= hxt=HFor("xt"), hyt=HFor("yt"); //获得函数xt和yt的句柄，保存在模块变量中
gleDrawScene[HFor("Scene")],stop(); //设置场景绘制函数后退出
Scene(::hxt,hyt,tmax,tmin)=
{
glClear[], //清除屏幕以及深度缓存
glLoadIdentity[], //重置视图
glTranslated[0,0,-20], //移动坐标，向屏幕里移动10个单元
glColor3d[0,0,1], //设置颜色
fgPlot[hxt,tmin,tmax,FG_Y,-1,2], //绘制一元函数图象
glColor3d[1,0,0], //设置颜色
fgPlot[hyt,tmin,tmax,FG_Y,-1,2] //绘制一元函数图象
};! l# {3 Y) t2 A4 N6 V$ E) l

作者: ExxNEN 时间: 2020-4-20 13:39
用1stOpt更简单：

Variable t=[0:0.2:30],x=1,y=0.5;
Plot x,y;
ODEFunction x'= -x^3-y, y'= x-y^3;& D1 _ Z6 g4 p

欢迎光临 EDA365电子论坛网 (https://bbs.eda365.com/)