如何從算法部分提取值？

Question

我正在使用拋物線PDE與代數方程混合，再加上所有這些方程。 我在OM（1.9.1）（「residualFcn [some number]」）中使用了Euler方法（Dassl太慢）和大容忍（用於快速模擬）和recive錯誤（對於這兩種類型）。問題是求解器可以解決非線性系統（數學上，系統是正確的）。 第一個問題是什麼類型的方法在OM中使用歐拉方法進行積分（顯式或隱式或Crank-Nicholson..or ...）？所以我試圖解決它的數值（顯式歐拉方法（在代碼下面的「新[N]」）（也許問題可以是CFL條件），但我有問題（樣本重建的具體採樣時間） 所以，第二個問題是指重現特定採樣時間的值？！ 在下面的代碼中有數組「a [3]」。 想法是針對每個「ts（採樣時間）」重建3個節點中的值。 這怎麼能 我怎樣才能從算法部分當前值（即節點）傳遞到公式部分？ 我有.txt文件中的所有值，我用Matlab來繪製它們，但我不知道如何通過它，即在「方程」部分或其他方式 同樣的問題是「new [N]」，針對特定的採樣時間，節點的繪圖函數（N）如果delta（t）/（delta（x））^ 2> = 0.5（delta（t）定義用戶，並且參考公式部分，則delta（x）如下面的代碼所示，並且空間離散化用於方程部分（經典前饋方法）），數值穩定性是否滿足？ 同樣的例子，但算法部分？ 問候

這裏是代碼：

model Euler1D 
import Modelica.Utilities.*; 
parameter Integer N=10; //50 
parameter Real Lp=1e-6; 
parameter Real deltax=1/(N-1)*Lp; 
Real a[3]; 
Real old[N]; 
Real new[N]; 
Real b; 
equation 
a[1]=if 
    (time>5) then 0 else time+5; 
a[2]=time; 
a[3]=2; 
when 
(sample(0,1)) then 
d=b; 
end when; 
algorithm 
// IN t=ts; 
when (sample(0,1)) then 
for i in 1:2 loop 
    b:=a[i]; 
Streams.print(String(time)+" "+String(a[i])+ "  "+String(b), "C:/Some_Path/text.txt"); 
end for; 
end when; 

// Another problem 
old[1]:=10; 
old [N]:=0; 
new[1]:=10; 
new [N]:=0; 
// Boundary 
for i in 2:N-1 loop 
old [i]:=10; 
new[i]:=10; 
end for; 

for dx in deltax:deltax:Lp-deltax loop // spatial discretization 

for i in 2:N-1 loop 
(new[i]):=(old[i]+0.5*(old[i + 1] +old[i-1]- 2*old[i])); 
//def:=def+abs(new[i]-old[i]); 
end for; 
for i in 2:N-1 loop 
old[i]:=new[i]; // switch the values 
end for; 

for i in 1:N loop 
Streams.print(String(time)+" "+ String(new[2]), "C:/Some_Path/Anel_Nodes.txt"); 
end for; 

annotation (uses(Modelica(version="3.2"))); 
end Euler1D; 

Answer 1

考慮您的意見，這是我會怎麼構建模型：

model Euler1D 
    parameter Integer N=10 "Spatial discretization"; 
    parameter Modelica.SIunits.Length L=1; 
protected 
    parameter Modelica.SIunits.Length dx=L/(N-1) "Segment size"; 
    Real c[N] "Solution variable c"; 
    Real J[N] "Solution variable J"; 
    Real dc_dx "Spatial derivative at x==L"; 
    Real d2c_dx2[N] "Second spatial derivative of c"; 
initial equation 
    der(c[2:N-1]) = zeros(N-2); 
equation 
    // Equations for spatial derivatives 
    d2c_dx2[2:N-1] = { (c[i+1]-2*c[i]+c[i-1])/(dx^2) for i in 2:N-1}; 
    dc_dx = (c[N]-c[N-1])/dx; 

    // Boundary conditions 
    c[1] = 0 "Dirichlet B.C."; 
    dc_dx = 1 "Neumann B.C."; 

    // PDE 
    J = c .* c "J = c^2"; 
    der(c) = d2c_dx2+J "PDE"; 
end Euler1D; 

初始方程避免大的瞬態在模擬的開始。如此有效地，我提出的模型將爲您提供穩定的狀態解決方案。但是你可以使事物變化（例如邊界條件）。

可能還存在一些錯誤，但我很難說，因爲我對這個系統沒有多少直覺，我也無法知道我使用的值是否是物理的有意義的。

但希望它概述瞭如何構建這樣一個模型。這個模型運行在Dymola。我沒有在OpenModelica中嘗試過。所有這些都說了，你總是用PDE來解決這個問題，即時間積分方案不能明確看到物理學產生的隱式約束（例如CFL條件）。通常，可變時間步長算法具有誤差估計，並且它們可以間接地看到由違反這些隱式約束而引入的不穩定性。但他們必須通過他們「摸索」。我不確定這將會成爲上述模型的問題。

如果你想一點點的興奮添加到您的解決方案，您可以更改公式

dc_dx = 1 "Neumann B.C."; 

...到...

dc_dx = 0.7+0.2*sin(time) "Neumann B.C."; 

...，你會得到一個不錯的時間變化解決我注意到，在玩這個遊戲的時候，如果你把梯度變得過於陡峭，模型會變得不穩定。無論這是基礎數學方程的真正解決方案還是由於上述問題導致的數值假象，我都不知道。所以如果你想要一個穩定的解決方案，你可以放入這個模型是有限的。