Matlab求解微分方程并画图

Question

Matlab小白，求大神帮忙写个程序，顺便学习一下，感恩。不出意外的话，图应该会长得像图2的样子。。。谢谢！

Answer 1

看标题以为你要求微分方程呐，结果是画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)

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