Какая разница в этом 2 кода?

Question

Какая разница в этом 2 кода?

Я написал два кода для решения дифференциальной задачи, оба с использованием подхода предиктора-корректора.

Первый предназначен для решения дифференциальной системы, так что в основном он использует вектор y с размером, равным номеру уравнения, которое нужно решить в системе, так что это в основном перезаписывается каждый шаг и не сохраняет шаг раньше.

Второй вместо этого использует вектор размером, равным количеству шагов, поэтому легко контролировать, какие элементы мы используем.

Я написал две разные (но в основном идентичные) функции: одна использовалась для первой процедуры, а вторая — для второй.

Я сделал много отладки, но я не понимаю, что я сделал не так.

Вот две разные версии решателя:

# include <iostream>
# include <functional>
# include <vector>
# include <string>
# include <cmath>
# include <fstream>
# include <valarray>using namespace std;std::vector<std::function<const double(const double t , std::valarray<double> u)>> func ;auto f2      = [](double t , double u) {return -20*u+20*sin(t)+cos(t) ; } ;

auto numFun  = [](double t , std::valarray<double> u) {return -20*u[0]+20*sin(t)+cos(t) ; } ;int main (){

double t0 = 0.0 ;
double tf = 2.5 ;
double dt = 0.01;const int N = (tf-t0)/dt;

func.push_back(numFun);

auto& f  =    func ;

double t1 ;
std::valarray<double> y1 ;
y1.resize(func.size()) ;   // size = 1
std::valarray<double> yp1 ;
yp1.resize(func.size()) ;  // size = 1std::vector<double> t2(N+1)  ;
std::vector<double> y2(N+1)  ;
std::vector<double> yp2(N+1) ;std::vector<double> u(1);
u.at(0) = 1.0 ;

t1  = t0 ;

for(auto i=0; i< u.size() ; i++)
{
y1[i]= u.at(i);
yp1[i] = u.at(i);
}ofstream fn("output1.dat", ios::out);for(auto i=1 ; i <=N ; i++ )
{

fn << t1 << " " ;
for(auto j =0 ; j < y1.size() ; j++ )
fn << y1[j] << " " ;
fn << endl ;

for(auto j=0 ; j < y1.size() ; j++)
{
yp1[j] = yp1[j] + dt * f[j](t1 , y1) ;

y1[j] += dt/2 * ( f[j](t1 , y1) + f[j]( t1+dt , yp1 ));

}
t1+=dt ;
}
fn.close();

/* --------------------          Vector version      --------------------------------- */

ofstream fn2("output2.dat", ios::out);

t2.at(0) = t0 ;

y2.at(0) = u.at(0) ;
yp2.at(0) = u.at(0) ;fn2 << t2.at(0) << " " ;
fn2 << y2.at(0) << " " ;
fn2 << endl ;

for(auto i=1; i <= N  ; i++ )
{
t2.at(i) = t2.at(i-1) + dt ;

// Predictor step (Exp Euler)

yp2.at(i) = y2.at(i-1) + dt * f2(t2.at(i-1), y2.at(i-1)) ;

y2.at(i) = y2.at(i-1) + dt/2 * ( f2(t2.at(i-1), y2.at(i-1)) + f2(t2.at(i), yp2.at(i)) ) ;fn2 << t2.at(i) << ' ' << y2.at(i) << " " ;

fn2 << endl;
}

fn2.close();return 0;
}

Я получил два разных результата в выходном файле:

-1

c++differential-equations

Решение

Другие решения

Чтобы первый метод работал для векторно-значных ODE, вам нужно разделить три шага на три цикла, соответствующие векторным операциям во втором методе. Обратите внимание, что второй метод содержит неявный промежуточный вектор в обновлении y1 так как правая часть сначала полностью вычисляется с использованием старого значения y1 перед копированием в y1,

  for(auto i=1 ; i <=N ; i++ )
{

fn << t1 << " " ;
for(auto j =0 ; j < y1.size() ; j++ )
fn << y1[j] << " " ;
fn << (exp(-20*t1) + sin(t1)) << " ";
fn << endl ;

for(auto j=0 ; j < y1.size() ; j++)
{
yp1[j] = y1[j] + dt * f[j](t1 , y1) ;

}
for(auto j=0 ; j < y1.size() ; j++)
{
yc1[j] = y1[j] + dt/2 * ( f[j](t1 , y1) + f[j]( t1+dt , yp1 ));
}
for(auto j=0 ; j < y1.size() ; j++)
{
y1[j] = yc1[j];
}
t1+=dt ;
}

3

Источник

Accepted Answer

попробуйте использовать 3 вектора таким образом!

         for(t=t0() ; t < tf() ; t+= dt()  )
{
f << t << ' ' ;
for(auto i=0 ; i< u.size() ; i++)
f << u[i] << " " ;
f << std::endl ;

for( auto i=0 ; i < u.size() ; i++)
{
up[i] = uc[i] + dt() * rhs.f[i](t , u) ;uc[i] += dt()/2 * ( rhs.f[i](t , u) + rhs.f[i]( t +dt() , up ) );

u[i] = uc[i] ;
}

}

Так что в основном переменная предиктора должна быть uc[i] + ... и не могло быть += как в стандартной формулировке Эйлера (в этом случае он написал yp1[j] знак равно yp1[j] + ... ) вместо yp1[j] = y[j] + ...

0