MATLAB：微分方程（组）数值解

news/2024/7/5 6:08:23

一、显式微分方程

 

clc,clear
tspan = [0:10]; 
y0 = 2;
[t1,y1] = ode23(@odefun_1,tspan,y0); %求数值解，精度相对低
[t2,y2] = ode113(@odefun_1,tspan,y0); %求数值解，精度相对高
yt = sqrt(tspan'+1)+1; %求精确解
subplot(1,2,1)
plot(t1,y1,'bo',t2,y2,'r*',tspan,yt,'g-')
legend('ode23求解','ode113求解','精确解','Location','best');
subplot(1,2,2)
y1t = abs(y1-yt); 
y2t = abs(y2-yt);
plot(t1,y1t,'bo',t2,y2t,'r*')
legend('ode23误差','ode113误差','Location','best');

function dydt = odefun_2(t,y)
    ft = 2*sin(t).*(t<4*pi) + 0.*(t>=4*pi); %定义分段函数f(t)
    gt = cos(t).*(t>=3.5*pi) + 0.*(t<3.5*pi); %定义分段函数g(t)
    dydt = [y(2) - ft;
           y(1)*gt-y(2)]; %微分方程组
end
ts = [0,20];
y0 = [1,2];
sol = ode23tb(@odefun_2,ts,y0)
subplot(1,2,1)
plot(sol.x,sol.y(1,:),'b-','LineWidth',1.5)
hold on
plot(sol.x,sol.y(2,:),'r--','LineWidth',1.5)
grid on
legend('{\ity}_1(t)','{\ity}_2(t)','Location','best')
legend('boxoff')
title('微分方程组解的曲线')
xlabel('t')
ylabel('y')
ysum = sum(sol.y);
subplot(1,2,2)
plot(sol.x,ysum,'b-','LineWidth',1.5)
hold on 
fplot(@(t)0*t,[0,20])
fh = @(t)sum(deval(sol,t));
[x01,y01] = fzero(fh,10)
[x02,y02] = fzero(fh,18)
plot(x01,y01,'r.','MarkerSize',15)
hold on
plot(x02,y02,'r.','MarkerSize',15)
grid on
legend('F(t) = {\ity}_1(t) + {\ity}_2(t)','y = 0的直线','Location','best');
legend('boxoff')

function dy = DyDtFun(t,y) 
    dy = zeros(3,1);
    dy(1) = y(2)*y(3);
    dy(2) = -y(1)*y(3);
    dy(3) = -0.51*y(1)*y(2);
end
[T,Y]=ode45(@DyDtFun,[0 12],[0 1 1]);
plot(T,Y(:,1),'b:',T,Y(:,2),'g--',T,Y(:,3),'r.-','LineWidth',1.5)
grid on
legend('y_1','y_2','y_3')
legend('boxoff')
title('显示微分方程组初值问题求解')

 

function dydt = odefun_3(t,y) %定义函数文件
    dydt = [y(2);
           (1-y(1)^2)*y(2)-y(1)];
end
[t,y]=ode45(@odefun_3,[0 20],[2 0]);
plot(t,y(:,1),'b--',t,y(:,2),'r:','LineWidth',1.5);
grid on
title('常微分方程的解');
xlabel('t') 
ylabel('y')
legend('y(t）','dy/dt','Location','best')
legend('boxoff')

 

function dxdt = apolloeq(t,x)
    mu = 1/82.45;
    mu1 = 1 - mu;
    r1 = sqrt((x(1)+mu).^2 + x(3).^2);
    r2 = sqrt((x(1)-mu1).^2 + x(3).^2);
    dxdt = [x(2);
            2*x(4) + x(1) - mu1*(x(1) + mu)/r1.^3 - mu*(x(1) - mu1)/r2.^3;
            x(4);
           -2*x(2) + x(3) - mu1*x(3)/r1.^3 - mu*x(3)/r2.^3];
end
x0 = [1.2;0;0;-1.04935751];
% [t,y] = ode45(@apolloeq,[0,20],x0);
[t,y] = ode45(@apolloeq,[0,20],x0,odeset('RelTol',1e-6));
plot(y(:,1), y(:,3),'LineWidth',2)
grid on
title('Apollo卫星的运动轨迹求解图像')

看课后案例分析。 

二、完全隐式微分方程

function dydt = odeifun_1(t,y,yp)
    dydt = t.*y.^2.*yp.^3 - y.^3.*yp.^2 + t.*(t.^2 + 1).*yp - t.^2.*y;
end
clc,clear
t0 = 1; 
y0 = sqrt(3/2); 
ypo = 0;
[y0,yp0] = decic(@odeifun_1,t0,y0,1,ypo,0)
%求得y0 = 1.2247 yp0 = 0.8165
[t,y] = ode15i(@odeifun_1,[1,10],y0,yp0,odeset('RelTol',1e-5));
plot(t,y)
hold on
plot(t,sqrt(t.^2+0.5),'r*') %绘制精确解散点
legend('y=ode15i','y=sqrt(t.^2+0.5)')

function dydt = odeifun_2(t,y,dy)
dydt = [dy(1)-y(2);
        dy(2)*sin(y(4))+dy(4)^2 + 2*y(1)*y(3)-y(1)*dy(2)*y(4);
        dy(3)-y(4);
        y(1)*dy(2)*dy(4)+cos(dy(4))-3*y(2)*y(3)];
end
t0 = 0; 
y0 = [1;0;0;1]; 
fix_y0 = ones(4,1); 
dy0 = [0 3 1 0]';
fix_dy0 = zeros(4,1); 
[y02,dy02] = decic(@odeifun_2,t0,y0,fix_y0,dy0,fix_dy0); 

[t y] = ode15i(@odeifun_2,[0 5],y02,dy02); 
plot(t,y(:,1),'k-',t,y(:,2),'r--', t,y(:,3),'g-.',t,y(:,4),'b:','linewidth',2)
grid on
legend('{\ity}_1(t)','{\ity}_2(t)','{\ity}_3(t)','{\ity}_4(t)','Location','best')
title('隐式微分方程组的解曲线')
xlabel('t')
ylabel('y')
legend('boxoff')

三、代数微分方程组（DAEs） 

ode15s用于求显式的微分方程组

y0 = [0.8;0.1;0.1]; %初值
M = [1 0 0;0 1 0;0 0 0]; %质量矩阵
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',[1e-6,1e-6,1e-6]);
tspan = [0 20];
[T,Y] = ode15s(@DAEsfun_1,tspan,y0,options);
plot(T,Y(:,1),'b-','linewidth',2)
hold on
plot(T,Y(:,2),'r-.','linewidth',2)
plot(T,Y(:,3),'g:','linewidth',2)
legend('{\ity}_1(t)','{\ity}_2(t)','{\ity}_3(t)','Location','best');
title('代数微分方程求解结果图')

ode15i用于求隐式的微分方程组 

function dydt = iDAEsfun_2(t,y,yp)
    % 注意与显式微分方程组定义的不同
    dydt = [yp(1)+0.3*y(1)+2*y(2).*sin(yp(3))+y(2).*y(4);
            yp(2)+y(2)+0.5*cos(yp(1)+y(3))+0.2*sin(0.6*t);
            yp(3)+0.2*y(1).*y(2)-exp(-yp(1));
            y(1)+y(2)-y(3)-y(4)-1];
end
y0 = [1;0.5;0.3;0.2]; % 初值
tspan = [0 5]; % 求值区间
fix_y0 = [0;1;1;1];
fix_yp0 = [0;0;0;0]; % 该组初值都可以改变，故全部为0
% fsolve
fh = @(y)[y(1)+0.4+sin(y(3)),y(2)+0.5+0.5*cos(y(1)+0.3),y(3)+0.1-exp(y(1))];
fsolve(fh,[1,1,1]) % -0.7596 -0.9481  0.3678
yp0 = [-0.8;-1;1;1];
% 首先求yp0，使得f(t,y,yp) = 0
[y02,yp02] = decic(@iDAEsfun_2,0,y0,fix_y0,yp0,fix_yp0)
[t,y] = ode15i(@iDAEsfun_2,tspan,y02,yp02);
plot(t,y(:,1),'k-',t,y(:,2),'r-.',t,y(:,3),'b:',t,y(:,4),'g--','linewidth',1.5)
grid on
legend('{\ity}_1(t)','{\ity}_2(t)','{\ity}_3(t)','{\ity}_4(t)','Location','best')
legend('boxoff')

四、延时/时滞微分方程 

 

显式微分方程 

function dydt = DDEsfun_1(t,y,Z)
    dydt = [0.5*Z(3,2) + 0.5*y(2).*cos(t);
            0.3*Z(1,1) + 0.7*y(3).*sin(t);
            y(2) + cos(2*t)];
end
lags = [1,3]; % 延迟常数向量
history = [0,0,1]; % 小于初值时的历史函数
tspan = [0,8];
sol = dde23(@DDEsfun_1,lags,history,tspan) % 调用dde23求解
ti = 0:0.01:8;
% yi = deval(sol,ti);
% plot(ti,yi(1,:),ti,yi(2,:),ti,yi(3,:))
plot(sol.x,sol.y(1,:),'r-','linewidth',2);
hold on
grid on
plot(sol.x,sol.y(2,:),'b-.','linewidth',2)
plot(sol.x,sol.y(3,:),'g-*','linewidth',1)
legend('{\ity}_1(t)','{\ity}_2(t)','{\ity}_3(t)','Location','best')
title('方程各解的曲线图')
legend('boxoff')

function dydt = DDEsfun_2(t,x,Z)
    dydt = [-2*x(2) - 3*Z(1,1);
            -0.05*x(1).*x(3) - 2*Z(2,2) + 2;
            0.3*x(1).*x(2).*x(3) + cos(x(1).*x(2)) + 2*sin(0.1*t.^2)];
end
lags = @(t,x)[t - 0.2*abs(sin(t)); t - 0.8];
history = zeros(3,1);
tspan = [0,10];
sol = ddesd(@DDEsfun_2,lags,history,tspan)
plot(sol.x,sol.y(1,:),'r--',sol.x,sol.y(2,:),'b-.',sol.x,sol.y(3,:),'g:','linewidth',2);
legend('{\ity}_1(t)','{\ity}_2(t)','{\ity}_3(t)','Location','best')
legend('boxoff')
grid on
title('变时间延迟微分方程各解的曲线图')

lags = @(t,x)[t - 0.2*abs(sin(t)); t - 0.8];
% history = zeros(3,1);
history = @(t,x)[sin(t+1); cos(t); exp(3*t)];
tspan = [0,10];
sol = ddesd(@DDEsfun_2,lags,history,tspan)
% ti = 0:0.01:10;
% yi = deval(sol,ti);
% plot(ti,yi(1,:),'r--',ti,yi(2,:),'b-.',ti,yi(3,:),'g:','linewidth',2)
plot(sol.x,sol.y(1,:),'r--',sol.x,sol.y(2,:),'b-.',sol.x,sol.y(3,:),'g:','linewidth',2)
legend('{\ity}_1(t)','{\ity}_2(t)','{\ity}_3(t)','Location','best')
legend('boxoff')
grid on
title('变时间延迟微分方程各解的曲线图')


function dxdt = DDEsfun_3(t,x,Z,ZP)
    A1 = [-13 3 -3; 106 -116 62; 207 -207 113];
    A2 = diag([0.02 0.03 0.04]);
    B = [0;1;2];
    U = 1;
    dxdt = A1*Z + A2*ZP + B*U;
end
history = zeros(3,1);
sol = ddensd(@DDEsfun_3,0.15,0.5,history,[0,15])
plot(sol.x,sol.y(1,:),'r--',sol.x,sol.y(2,:),'b-.',sol.x,sol.y(3,:),'g:','linewidth',2)
legend('{\ity}_1(t)','{\ity}_2(t)','{\ity}_3(t)','Location','best')
legend('boxoff')
grid on
title('中立型延迟微分方程组解曲线')

五、微分方程边值问题

 solinit: solve initial.

 

function dydt = bvpfun_1(t,y)
    dydt = [y(2);
            2*y(2).*cos(t) - y(1).*sin(4*t) - cos(3*t)];
end
function res = bcfun_1(ya,yb)
    res = [ya(1) - 1;
           yb(1) - 2];
end
function yint = initfun_1(t)
    %  y(0) = 1, y(4) = 2 --> (0,1),(4,2)
    yint = [1/4*t + 1;
            1/4]
end
function dfdy = jacb(t,y)
    dfdy = [0, 1;
            -sin(4*t), 2*cos(t)];
end
solinit = bvpinit(linspace(0,4,10),@initfun_1)
opts = bvpset('FJacobian',@jacb,'RelTol',1e-4,'AbsTol',1e-8,'Stats','on');
sol = bvp4c(@bvpfun_1,@bcfun_1,solinit,opts)
sol5 = bvp5c(@bvpfun_1,@bcfun_1,solinit,opts)
% bvp5c精度比bvp4c高
ti = 0:0.01:4;
yi = deval(sol,ti);
plot(ti,yi(1,:),ti,yi(2,:))
grid on
legend('{\ity}_1(t)','{\ity}_2(t)','Location','best')
legend('boxoff')
title('常微分方程边界问题解的曲线图')


http://lihuaxi.xjx100.cn/news/2125457.html

相关文章

linux 使用kubeadm安装部署K8S

linux 使用kubeadm安装部署K8S

上一节&#xff0c;安装了docker-20.10.14版本 安装kubeadm&#xff0c;kubelet和kubectl 由于版本更新频繁&#xff0c;这里指定版本号部署 yum install kubelet-1.23.4-0 kubeadm-1.23.4-0 kubectl-1.23.4-0 -y Loaded plugins: fastestmirror Repository epel is listed m…
阅读更多...
【MySQL】聊聊自增id用完怎么办？

【MySQL】聊聊自增id用完怎么办？

在实际的开发中&#xff0c;一般都会将数据存储到数据库中&#xff0c;在设计表的时候&#xff0c;其实id如果达到最大值的话&#xff0c;会出现什么问题。其实主要分两种情况&#xff0c;一种是设置了主键id&#xff0c;另一种没有设置主键id。 表定义自增值id create table…
阅读更多...
【蓝桥杯】tarjan算法

【蓝桥杯】tarjan算法

一.概述 Tarjan 算法是基于DFS的算法&#xff0c;用于求解图的连通性问题。 Tarjan 算法可以在线性时间内求出&#xff1a; 无向图&#xff1a; 割点与桥双连通分量 有向图&#xff1a; 强连通分量必经点与必经边 1.割点&#xff1a; 若从图中删除节点 x 以及所有与 x 关联的…
阅读更多...
【机器学习】原理+实例， 一文掌握 精确率、召回率与F1分数，再也不迷路。

【机器学习】原理+实例， 一文掌握 精确率、召回率与F1分数，再也不迷路。

精确率、召回率与F1分数 1、引言2、定义2.1 精确率2.2 召回率2.3 F1分数2.4 二分类任务2.5 代码示例 3、总结 1、引言 小屌丝&#xff1a;鱼哥&#xff0c;你在给我详细的唠叨唠叨 精确率&#xff0c;召回率和F1分数。 小鱼&#xff1a;这…这不是机器学习基本知识吗 小屌丝&a…
阅读更多...
Tomcat 下载以及安装

Tomcat 下载以及安装

Tomcat安装及配置教程主要分为四步&#xff1a; 步骤一&#xff1a;首先确认自己是否已经安装JDK 1. cmd&#xff1a;查看java的版本 步骤二&#xff1a;下载安装Tomcat 1. 下载tomcat :Apache Tomcat - Welcome! 2. 选择对应的tomcat版本&#xff1a; 3. 进行安装&#…
阅读更多...
【数据分享】1929-2023年全球站点的逐日平均压力（Shp\Excel\免费获取）

【数据分享】1929-2023年全球站点的逐日平均压力（Shp\Excel\免费获取）

气象数据是在各项研究中都经常使用的数据&#xff0c;气象指标包括气温、风速、降水、能见度等指标&#xff0c;说到气象数据&#xff0c;最详细的气象数据是具体到气象监测站点的数据&#xff01; 有关气象指标的监测站点数据&#xff0c;之前我们分享过1929-2023年全球气象站…
阅读更多...
前后端实时数据通信

前后端实时数据通信

实现前后端实时数据转换通常涉及到以下几个步骤&#xff1a; 后端提供数据转换接口。 前端实时数据获取。 前端实时数据转换。 前端实时展示转换后数据。 以下是一个简单的例子&#xff0c;假设后端提供了一个接口来转换某种数据格式&#xff0c;前端使用JavaScript和WebS…
阅读更多...
Redis缓存预热，该如何实现

Redis缓存预热，该如何实现

一、什么是缓存预热 缓存预热是一种在程序启动或缓存失效之后&#xff0c;主动将热点数据加载到缓存中的策略。这样&#xff0c;在实际请求到达程序时&#xff0c;热点数据已经存在于缓存中&#xff0c;从而减少了缓存穿透和缓存击穿的情况&#xff0c;也缓解了SQL服务器的压力…
阅读更多...
最新文章

Copyright @ 2022~2023