刚性方程的隐式欧拉

Question

问题：求解刚性微分方程
方法：隐式欧拉
计划：我通过求解非线性方程使用正割法计算下一个'y'。我的函数是 dy/dx = sin(x+y)

有正确的解决方案。我用的是牛顿法

main.m

h=0.01;
x(1)=0;
y_expl(1)=0;  
y_impl(1)=0+h;  
dy(1)=0;
eps=1.0e-6;

for i=1:1000

    x(i+1)=x(i)+h;

    y_impl(i+1)=newton(x(i),y_impl(i),y_impl(i));

    y_expl(i+1)=y_expl(i)+h*f(x(i),y_expl(i));

end

plot(x,y_impl,'r',x,y_expl,'b')
legend('Implicit Euler','Explicit Euler');

newton.m

function [ yn ] = newton( x,y,yi )

    eps=1.0e-6;
    err=1;
    step=0;
    step_max=100;
    h=0.01;
    xn=x+h;

    while (err > eps) && (step < step_max)

        step=step+1;

        yn=y-(F(xn,y,yi,h))/(J(xn,y,h));

        err=abs(y-yn)/(abs(yn)+1.0e-10);

        y=yn;

    end

end

f.m

function [ res ] = f( x,y )

res = sin(x+y);

end

G.m

function [ res ] = J( xn,y,h )

res = h*f(xn,y)-1;

end

F.m

function [ res ] = F( a,y,yn,h )

res = h*f(a,y)-y+yn;

end

感谢关注

Answer 1

问题是您不应该求解 F(x,y)=0，而是求解由隐式欧拉步骤 y=y0+h*F(x,y) 得出的方程。从而定义

function [res] = G(x,y,y0,h)
    res = y - y0 - h*F(x,y)
end

并使用牛顿法或正割法 G。

一般代码评论：

matlab中没有x(0)
隐式欧拉是一种单步法，不需要为索引 2 和 3 进行初始化。
x 值的迭代是 x(i+1)=x(i)+h

在正割法中：已知值为 x0=x(i)、y0=y(i) 和 h。

开始循环之前需要 x1=x0+h 和 y1 的初始值。这可以作为显式欧拉步骤的结果，y1=y0+h*F(x0,y0) 作为预测变量。正割法用作校正器。
如果G的值是分开计算的，那么代码的可读性更好。注意G(x1,y,y0,h)中的变量是y，其他都是固定参数。因此计算 G0=G(x1,y0,y0,h) 和 G1=G(x1,y1,y0,h) 用于割线公式 y2=y1-G1*(y1-y0)/(G1-G0) 或更对称的 y2=(y0*G1-y1*G0)/(G1-G0).
原则上，您可以通过调用
使用具有接口 secant(func, a, b, tol) 的通用正割方法
```
x(i+1) = x(i)+h;
y(i+1) = secant(@(y) G(x(i+1),y,y(i),h), y(i), y(i)+h*F(x(i),y(i)), delt)
```

刚性方程的隐式欧拉

Implicit Euler for stiff equation

matlab

numerical-methods

numerical-integration

differential-equations