计算科学 - ode45 在这种情况下的用法？ - 吾爱随笔录

ode45 在这种情况下的用法？

计算科学 matlab 颂

2021-11-29 05:18:09

背景是我正在解决数值分析中的一个问题，我似乎已经解决了非线性系统的一半，现在到了我应该应用求解微分方程的部分。

https://math.stackexchange.com/questions/467164/how-to-proceed-solving-this-problem

到目前为止的代码是

z=[0:500:4000 5000:1000:12000];
data=[5050 4980 4930 4890 4870 4865 4860 4860 4865 4875 4885 4905 4920 4935 4950 4970 4990];
fun=@(p)(4800 + p(1))*ones(size(z)) +p(2)/1000*z+p(3)*exp(p(4)/1000*z)-data;
x0=[0 0 0 -1];
opt = optimset('MaxFunEvals',1000);
p=lsqnonlin(fun,x0,[],[],opt);
fitf=@(t)(4800 + p(1))*ones(size(t)) + p(2)/1000*t+ p(3)*exp(p(4)/1000*t);
tt=linspace(0,12000,1000);
plot(z,data,'r-',tt,fitf(tt),'b-');

方程组在数学问题中说明。现在我想我应该定义一个类似的函数

function f=fsys(x, z)
    f= ...

但首先我想我必须将系统从二阶重写为一阶，这只是数学理论，对吧？所以我首先将系统重写为一阶，然后在 matlab 中定义函数，最后类似于

xend=...; z0=...
[X,Z]=ode45(@fsys,[0, xend],z0); zs=Z(end), plot(X,Z,'--')

你能帮我定义 fsys，告诉我如何继续，确认或拒绝我走在正确的路上吗？如果您不想帮助我，但提示如何学习 matlab 来解决这个“非常简单的数学”，有什么提示吗？

我想我可以用通常的方式改造系统：

u = z'
u' = z''

系统变成

u'=-q0*c'(z)/((c(z)^3)

由于 z 是 x 的函数，u 和 u' 也是 x 的函数，因此该系统直到成为一阶系统，但是如何在 matlab 中编程解决方案？如何使用 u=tan(b0)？

1个回答

我取并得到其他系数为，，。所以和将是 $p4=1$ $p_{1}=5.66$ $p_{2}=14.58$ $p_{3}=257.4$ $c(z)$ $c'(z)$

$c(z)=4800+5.66+(14.58)\frac{z}{1000}+(257.4)exp(-\frac{z}{1000})$ 和

$c'(z)=(1287exp(\frac{z}{1000}))/5000+729/50000$

让和比我们可以写以下微分方程。 $z'=u$ $z''=u'$

\frac{d^{2} z}{d x^{2}} = {(\frac{c (2000)}{c o s (7.8)})}^{2} \frac{c^{'} (z)}{c (z)^{3}}, z (0) = 2000, z^{'} (0) = t a n (7.8)

$\frac{d^{2}z}{dx^{2}}=\left(\frac{c(2000)}{cos(7.8)}\right)^{2}\frac{c'(z)}{c(z)^{3}}, \qquad z(0)=2000,\qquad z'(0)=tan(7.8)$

作为

\begin{aligned} z^{'} & = u, z (0) = 2000 \\ u^{'} & = {(\frac{c (2000)}{c o s (7.8)})}^{2} \frac{c^{'} (z)}{c (z)^{3}}, u (0) = t a n (7.8) \end{aligned}

$\begin{align*} z' & =u,\qquad z(0)=2000\\ u' & =\left(\frac{c(2000)}{cos(7.8)}\right)^{2}\frac{c'(z)}{c(z)^{3}},\qquad u(0)=tan(7.8) \end{align*}$

所以问题的代码应该是

    function dZ=sys(x,Z)
    c=@(z)4800 + 5.66 + (14.58)*z/1000+ (257.4)*exp(-z/1000); % c(z)
    c=c(2000);

% Z(1):=z
% Z(2):=u
    dZ=zeros(2,1);    % a column vector
    dZ(1)=Z(2);
    dZ(2)=(c/cosd(7.8))^2*((1287*exp(Z(1)/1000))/5000 + 729/50000)/...
        (4800 + 5.66 + (14.58)*Z(1)/1000+ (257.4)*exp(-Z(1)/1000))^3;
    end

我们可以通过以下代码调用上述函数并绘制解决方案

>> x=0:0.5:3000;
>> [X,Z]=ode45(@sys,x,[2000 tand(7.8)]);
>> plot(X,Z(:,1),'r',X,Z(:,2),'b')  %Z(:,1) is z(x) and Z(:,2) is z'(x).

你只看到下面的 z 在此处输入图像描述

其它你可能感兴趣的问题

上一篇有限差分计算中的误差范数下一篇在并行 PETSC 程序中创建局部向量