deval函数_Matlab求解微分方程并画图

展开全部

^

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

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