Задача (на Коши) в Mathlab

[Гост]

Здравейте, имам проблем с тези две задачи в Mathlab.
Всяко мнение ще ми бъде от полза!

[Anubis]

Създаваме m-файл с име ode22, който:
1. намира точното решение чрез функцията dsolve();
2. пресмята приближеното решение по явния метод на Ойлер с грешка от първи ред; това означава, че разликата между точното и приближеното решение е от порядъка на стъпката h;
3. пресмята приближеното решение по метода на Кранк-Никълсън с грешка от втори ред, което показва, че разликата между точното и приближеното решение е от порядъка на квадрата на избраната стъпка h; тъй като сме избрали h=0.01, разликата между решенията няма да надминава 0.0001.

Код: Избери целия код

clear all
clc
tic
%===================================
% Solving the equation exactly
%===================================
sol = ...
     dsolve('Dy=-x^3-x*y+4,y(0)=-3','x');
 simplify(sol);
 display('Exact solution: ')
 pretty(sol)
%===================================
% Grid
%===================================
x0 = 0; xend = 2;
%===================================
% Gird for the first numerical solution
%===================================
h1 = 0.5;
x1 = x0 : h1 : xend;
N1 = length(x1);
y1 = zeros(1,N1);
%===================================
% Initial condition
%===================================
y1(1) = -3;
%===================================
% Explicit Euler method with error O(h)
%===================================
for i = 1 : N1-1
    y1(i+1) = y1(i) + h1 .* ...
        (-(x1(i))^3 - x1(i) .* y1(i) + 4);
end
%===================================
% Grid for the second numerical solution
%===================================
h2 = 0.4;
x2 = x0 : h2 : xend;
N2 = length(x2);
y2 = zeros(1,N2);
%===================================
% Initial condition
%===================================
y2(1) = -3;
%===================================
% Explicit Euler method with error O(h)
%===================================
for i = 1 : N2-1
    y2(i+1) = y2(i) + h2 .* ...
        (-(x2(i))^3 - x2(i) .* y2(i) + 4);
end
%===================================
% Grid for the third numerical solution
%===================================
h3 = 0.1;
x3 = x0 : h3 : xend;
N3 = length(x3);
y3 = zeros(1,N3);
%===================================
% Initial condition
%===================================
y3(1) = -3;
%===================================
% Explicit Euler method with error O(h)
%===================================
for i = 1 : N3-1
    y3(i+1) = y3(i) + h3 .* ...
        (-(x3(i))^3 - x3(i) .* y3(i) + 4);
end
%===================================
% Finer grid for the exact solution
%===================================
xi = x0 : 0.01 : xend;
yi = subs(sol,'x',xi);
%===================================
% Using implicit Euler method (Crank-Nicolson)
% with error O(h^2)
%===================================
h = 0.01;
x = x0 : h : xend;
N = length(x);
sOld = zeros(1,N);
sNew = zeros(1,N);
sNew(1) = -3;
for i = 1 : N-1
    sNew(i+1) = sNew(i) + ...
        (h/2) .* ...
        ((-(x(i))^3-x(i)*sNew(i)+4) + ...
        (-(x(i+1))^3-x(i+1)*sOld(i+1)+4));
end
eps = 1.0e-06;
s = 1;
while(abs(max(sNew-sOld))>eps)
    sOld = sNew;
    sNew(1) = -3;
    for i = 1 : N-1
        sNew(i+1) = sNew(i) + ...
            (h/2) .* ...
            ((-(x(i))^3-x(i)*sNew(i)+4) + ...
            (-(x(i+1))^3-x(i+1)*sOld(i+1)+4));
    end
    s = s+1;
end
y = sNew;
%===================================
% Plot the solutions
%===================================
figure(1)
plot(x1,y1,'b--','LineWidth',3)
hold on
grid on
plot(x2,y2,'g--','LineWidth',3)
plot(x3,y3,'m--','LineWidth',3)
plot(xi,yi,'r','LineWidth',3)
plot(x,y,'cyan--','LineWidth',3)
set(gca,'FontName','Times','FontSize',13)
xlabel('$$ x \in [0; \, 2] $$','interpreter','latex')
ylabel('\it{Solution}')
legend('\it{y_{1}}','\it{y_{2}}','\it{y_{3}}','\it{Exact}','\it{Crank-Nicolson}')
toc


Трябва да се отбележи, че в практиката много по-често използваните методи за приближено решаване на диференциални уравнения са неявните. За изчисляване на координатите на приближеното решение се прилага итерационна процедура, както се вижда от скрипта ми. Необходими са само няколко итерации (s=7) за достигане на точност 0.000001 между две последователни приближения.