学习笔记-Matlab算法篇-方程求解_微分方差数学建模-CSDN博客

本文链接：https://blog.csdn.net/seek97/article/details/108338145

方程求解

01微分方差建模

介绍：微分方程建模是数学建模的重要方法，因为许多实际问题的数学描述将导致求解微分方程的定解问题。把形形色色的实际问题化成微分方程的定解问题，大体上可以按以下几步:

1.根据实际要求确定要研究的量(自变量、未知函数、必要的参数等)并确定坐标系。

2.找出这些量所满足的基本规律(物理的、几何的、化学的或生物学的等等)。

3.运用这些规律列出方程和定解条件。

列方程常见方法有：

1.按规律直接列方程

2.微元分析法与任意区域上取积分的方法

3.模拟近似法

例子(e01)：人口增长模型。据考古学家论证，地球上出现生命距今已有 20 亿年，而人类的出现距今却不足 200万年。纵观人类人口总数的增长情况，我们发现：1000 年前人口总数为 2.75 亿。经过漫长的过程到 1830 年，人口总数达 10 亿，又经过 100 年，在 1930 年，人口总数达 20亿；30 年之后，在 1960 年，人口总数为 30 亿；又经过 15 年，1975 年的人口总数是40 亿，12 年之后即 1987 年，人口已达 50 亿。

我们自然会产生这样一个问题：人类人口增长的规律是什么？如何在数学上描述这一规律。

Malthus模型：1789 年，英国神父Malthus在分析了一百多年人口统计资料之后，提出了 Malthus模型。

模型假设：

(1) 设𝑥(𝑡)x(t)表示t时刻的人口数，且𝑥(𝑡)x(t)连续可微。

(2) 人口的增长率r是常数（增长率=出生率-死亡率）

(3) 人口的数量变化是封闭的。

t=1961:1970;
y=3.06*10^9*exp(t-1961);
figure,
subplot(2,1,1)
plot(t,y,'b--*','linewidth',5);
xlabel('年份/年');
ylabel('人口/口');

subplot(2,1,2)
t=1961:2000;
y=3.06*10^9*exp(t-1961);
plot(t,y,'r-+','linewidth',5);
xlabel('年份/年');
ylabel('人口/口');

用matlab求解微分方程：

clear global;clc
syms r t0 xm x0;
xt=dsolve('Dx=r*(1-x/xm)*x','x(t0)=x0','t');
disp(xt);

%%对x做2阶导数：D2x

结果：

-xm/(exp(xm*((log(1 - xm/x0) + r*t0)/xm - (r*t)/xm)) - 1)

https://www.un.org/zh/sections/issues-depth/population/

02常微分方差解法

Matlab求解：Matlab的工具箱提供了几个解非刚性常微分方程的功能函数，如 ode45，ode23，ode113，其中ode45采用四五阶RK方法，是解非刚性常微分方程的首选方法，ode23采用二三阶RK方法，ode113采用的是多步法，效率一般比ode45高。

Matlab实现改进的欧拉方程：

function [x,y]=eulerpro(fun,x0,xfinal,y0,n)
    if nargin<5
        n=50;
    end
    h=(xfinal-x0)/n;
    x(1)=x0;y(1)=y0;
    for i=1:n
        x(i+1)=x(i)+h;
        y1=y(i)+h*feval(fun,x(i),y(i));
        y2=y(i)+h*feval(fun,x(i+1),y1);
        y(i+1)=(y1+y2)/2;
    end
end

function f=doty(x,y)
    f=-2*y+2*x.^2+2*x;
end

%% 欧拉方程求解
[x,y]=eulerpro('doty',0,0.5,1);
plot(x,y);

结果：

% RK方法求解
[x,y]=ode45('doty',[0,0.5],1);
figure,
plot(x,y,'r--*');title('ODE45')
figure,
[x,y]=ode23('doty',[0,0.5],1);
plot(x,y,'b-+');title('ODE23')
figure,
[x,y]=ode113('doty',[0,0.5],1);
plot(x,y,'g-.o');title('ODE113')

结果：

%% 求解对比
[x,y]=ode45('doty',[0,0.5],1);
figure
plot(x,y,'r--*');
hold on
[x,y]=ode23('doty',[0,0.5],1);
plot(x,y,'b-+');
[x,y]=ode113('doty',[0,0.5],1);
plot(x,y,'g-.o');

[x,y]=eulerpro('doty',0,0.5,1);
plot(x,y,'-->');
legend('ODE45','ODE23','ODE113','EULER-PRO');
grid on
xlabel('X')
ylabel('Y')

结果：

刚性常微分方程：Matlab的工具箱提供了几个解刚性常微分方程的功能函数，如ode15s，ode23s，ode23t，ode23tb，这些函数的使用同非刚性微分方程的功能函数。

常微分方程解析解：dsolve函数

%% 常微分防蹭解析解
y=dsolve('x^2+y+(x-2*y)*Dy=0','x');
disp(y);

结果：

>> e04
 x/2 + ((4*x^3)/3 + x^2 + C1)^(1/2)/2
 x/2 - ((4*x^3)/3 + x^2 + C1)^(1/2)/2

%% 带有初值条件的常微分方程解析解
y=dsolve('D3y-D2y=x','y(1)=8','Dy(1)=7','D2y(2)=4','x');
disp(y);

结果：

>> e05
x*((exp(-1)*(19*exp(1) - 14))/2 - 1) + 7*exp(-2)*exp(x) - x^2/2 - x^3/6 + (exp(-1)*(19*exp(1) - 14))/2 - (exp(-1)*(25*exp(1) - 21))/3 - 1

%% 求解常微分方程组
clc,clear
equ1='D2f+3*g=sin(x)';
equ2='Dg+Df=cos(x)';
[general_f,general_g]=dsolve(equ1,equ2,'x')
[f,g]=dsolve(equ1,equ2,'Df(2)=0,f(3)=3,g(5)=1','x')

结果：

general_f =
 
C5 + sin(x) + (3^(1/2)*exp(3^(1/2)*x)*(C6 + (exp(-3^(1/2)*x)*(cos(x) - 3^(1/2)*sin(x)))/4))/3 - (3^(1/2)*exp(-3^(1/2)*x)*(C7 + (exp(3^(1/2)*x)*(cos(x) + 3^(1/2)*sin(x)))/4))/3
 
 
general_g =
 
(3^(1/2)*exp(-3^(1/2)*x)*(C7 + (exp(3^(1/2)*x)*(cos(x) + 3^(1/2)*sin(x)))/4))/3 - (3^(1/2)*exp(3^(1/2)*x)*(C6 + (exp(-3^(1/2)*x)*(cos(x) - 3^(1/2)*sin(x)))/4))/3
 
 
f =
 
sin(x)/2 - sin(3)/(2*(exp(4*3^(1/2)) - exp(2*3^(1/2)) + 1)) + 3/(exp(4*3^(1/2)) - exp(2*3^(1/2)) + 1) - (2*exp(2*3^(1/2)))/(exp(4*3^(1/2)) - exp(2*3^(1/2)) + 1) + (3*exp(4*3^(1/2)))/(exp(4*3^(1/2)) - exp(2*3^(1/2)) + 1) - (exp(5*3^(1/2))*exp(-3^(1/2)*x))/(exp(6*3^(1/2)) + 1) + (exp(2*3^(1/2))*sin(3))/(2*(exp(4*3^(1/2)) - exp(2*3^(1/2)) + 1)) - (exp(2*3^(1/2))*sin(5))/(2*(exp(4*3^(1/2)) - exp(2*3^(1/2)) + 1)) - (exp(4*3^(1/2))*sin(3))/(2*(exp(4*3^(1/2)) - exp(2*3^(1/2)) + 1)) - (exp(3^(1/2)*x)*exp(3^(1/2)))/(exp(6*3^(1/2)) + 1) + (exp(3^(1/2)*x)*exp(3^(1/2))*sin(5))/(2*(exp(6*3^(1/2)) + 1)) + (3^(1/2)*cos(2)*exp(3^(1/2)))/(6*(exp(4*3^(1/2)) - exp(2*3^(1/2)) + 1)) + (exp(5*3^(1/2))*exp(-3^(1/2)*x)*sin(5))/(2*(exp(6*3^(1/2)) + 1)) - (3^(1/2)*cos(2)*exp(3*3^(1/2)))/(6*(exp(4*3^(1/2)) - exp(2*3^(1/2)) + 1)) - (3^(1/2)*cos(2)*exp(-2*3^(1/2))*exp(3^(1/2)*x))/(6*(exp(6*3^(1/2)) + 1)) + (3^(1/2)*cos(2)*exp(8*3^(1/2))*exp(-3^(1/2)*x))/(6*(exp(6*3^(1/2)) + 1))
 
 
g =
 
sin(x)/2 + (exp(5*3^(1/2))*exp(-3^(1/2)*x))/(exp(6*3^(1/2)) + 1) + (exp(3^(1/2)*x)*exp(3^(1/2)))/(exp(6*3^(1/2)) + 1) - (exp(3^(1/2)*x)*exp(3^(1/2))*sin(5))/(2*(exp(6*3^(1/2)) + 1)) - (exp(5*3^(1/2))*exp(-3^(1/2)*x)*sin(5))/(2*(exp(6*3^(1/2)) + 1)) + (3^(1/2)*cos(2)*exp(-2*3^(1/2))*exp(3^(1/2)*x))/(6*(exp(6*3^(1/2)) + 1)) - (3^(1/2)*cos(2)*exp(8*3^(1/2))*exp(-3^(1/2)*x))/(6*(exp(6*3^(1/2)) + 1))

02偏微分方差解法

将只含有未知多元函数及其偏导数的方程，称之为偏微分方程。方程中出现的未知函数偏导数的最高阶数称为偏微分方程的阶。

如果方程中对于未知函数和它的所有偏导数都是线性的，这样的方程称为线性偏微分方程，否则称它为非线性偏微分方程。

初始条件和边界条件称为定解条件，未附加定解条件的偏微分方程称为泛定方程。

对于一个具体的问题，定解条件与泛定方程总是同时提出。定解条件与泛定方程作为一个整体，称为定解问题。

偏微分方程的差分解法：差分方法又称为有限差分方法或网格法，是求偏微分方程定解问题的数值解中应用最广泛的方法之一。它的基本思想是：先对求解区域作网格剖分，将自变量的连续变化区域用有限离散点（网格点）集代替；将问题中出现的连续变量的函数用定义在网格点上离散变量的函数代替；通过用网格点上函数的差商代替导数，将含连续变量的偏微分方程定解问题化成只含有限个未知数的代数方程组（称为差分格式）。如果差分格式有解，且当网格无限变小时其解收敛于原微分方程定解问题的解，则差分格式的解就作为原问题的近似解（数值解）。因此，用差分方法求偏微分方程定解问题一般需要解决以下问题：

（i）选取网格；

（ii）对微分方程及定解条件选择差分近似，列出差分格式；

（iii）求解差分格式；

（iv）讨论差分格式解对于微分方程解的收敛性及误差估计。

function e07
%************************************
%求解一维热传导偏微分方程的一个综合函数程序
%************************************
m=0;
x=linspace(0,1,20); %xmesh
t=linspace(0,2,20); %tspan
%************
%以 pde 求解
%************
sol=pdepe(m,@ex20_1pdefun,@ex20_1ic,@ex20_1bc,x,t);
u=sol(:,:,1); %取出答案
%************
%绘图输出
%************
figure(1)
surf(x,t,u)
title('pde 数值解')
xlabel('位置 x')
ylabel('时间 t' )
zlabel('数值解 u')
%*************
%与解析解做比较
%*************
figure(2)
surf(x,t,exp(-t)'*sin(pi*x));
title('解析解')
xlabel('位置 x')
ylabel('时间 t' )
zlabel('数值解 u')
%*****************
%t=tf=2 时各位置之解
%*****************
figure(3)
M=length(t); %取终点时间的下表
xout=linspace(0,1,100); %输出点位置
[uout,dudx]=pdeval(m,x,u(M,:),xout);
plot(xout,uout); %绘图
title('时间为 2 时,各位置下的解')
xlabel('x')
ylabel('u')
%******************
%pde 函数
%******************
function [c,f,s]=ex20_1pdefun(x,t,u,dudx)
c=pi^2;
f=dudx;
s=0;
%******************
%初始条件函数
%******************
function u0=ex20_1ic(x)
u0=sin(pi*x);
%******************
%边界条件函数
%******************
function [pl,ql,pr,qr]=ex20_1bc(xl,ul,xr,ur,t)
pl=ul;
ql=0;
pr=pi*exp(-t);
qr=1;

结果：

练习(p01)：求解以下联立的偏微分方程：

function p01
%***************************************
%求解一维偏微分方程组的一个综合函数程序
%***************************************
m=0;
x=[0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1];
t=[0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];
%*************************************
%利用 pdepe 求解
%*************************************
sol=pdepe(m,@ex20_2pdefun,@ex20_2ic,@ex20_2bc,x,t);
u1=sol(:,:,1); %第一个状态之数值解输出
u2=sol(:,:,2); %第二个状态之数值解输出
%*************************************
%绘图输出
%*************************************
figure(1)
surf(x,t,u1)
title('u1 之数值解')
xlabel('x')
ylabel('t')
%
figure(2)
surf(x,t,u2)
title('u2 之数值解')
xlabel('x')
ylabel('t')
%***************************************
%pde 函数
%***************************************
function [c,f,s]=ex20_2pdefun(x,t,u,dudx)
c=[1 1]';
f=[0.024 0.170]'.*dudx;
y=u(1)-u(2);
F=exp(5.73*y)-exp(-11.47*y);
s=[-F F]';
%****************************************
%初始条件函数
%****************************************
function u0=ex20_2ic(x)
u0=[1 0]';
%****************************************
%边界条件函数
%****************************************
function [pl,ql,pr,qr]=ex20_2bc(xl,ul,xr,ur,t)
pl=[0 ul(2)]';
ql=[1 0]';
pr=[ur(1)-1 0]';
qr=[0 1]';

结果：