Тема 3. Интерполяция

Покажем, как решить задачу интерполяции на Scilab. Используются приведенные ниже исходные данные.

xbasc()

clc

function y=G(c,z)

y=z(2)-c(1)-c(2)*z(1)-c(3)*z(1)^2-c(4)*z(1)^3

endfunction

//Исходные данные

x=[-2 0 2 4];y=[-12 -10 -8 42];

//Формирование матрицы исходных данных

z=[x;y];//та же буква, что и в функции y=G(c,z)

//Вектор начальных приближений

c=[0;0;0;0];//нулей столько, сколько искомых коэффициентов. Это нач. приближение

//Решение задачи

[a,err]=datafit(G,z,c)

//Построение графика экспериментальных данных

plot2d(x,y,-4);

//Построение графика подобранной функции

t=-2:.1:4;

Ptc=a(1)+a(2)*t+a(3)*t^2+a(4)*t^3;

plot2d(t,Ptc);

t=1;

Ptc=a(1)+a(2)*t+a(3)*t^2+a(4)*t^3

err = 5.191D-12

a =

- 10.000001

- 2.9999989

0.0000006

0.9999998

Ptc =-12.

Здесь вектор а – вектор коэффициентов интерполяционного полинома, Ptc - решение задачи интерполяции, err - сумма площадей квадратов отклонений (то есть 0 в нашем случае).

График интерполяционного полинома

Второй вариант

clc

xbasc()

x=[6 8 10 12 14];y=[2.5 4.9 8 12.1 16.9 ];nst=51;

n=length(x);a=[];b=[];c=[];

for i=1:n

for j=1:n

a(i,j)=sum(x.^(i+j-2));

end

b(i)=sum(x.^(i-1).*y);

end

c=inv(a)*b

function z=f(t)

z=0;

for i=1:n

z=z+t.^(i-1).*c(i);

end

endfunction

plot2d(x,y,-4)

t=min(x)-.1:.01:max(x)+.1;plot2d(t,f(t))

r=min(x)+nst/30,z1=f(r),plot(r,z1,"*")

err=sum((y-f(x)).^2)

c = - 12.5 5.425 - 0.80625 0.0625 - 0.0015625

r = 7.7 z1 = 4.5105873 err = 6.351D-16

Численное интегрирование

Пример.

Вычислить

-->integrate('x/(x+1)','x',0,4) // ( ans = 2.3905621 )

Тема 4. Численные методы решения обыкновенных дифференциальных уравнений Решение дифференциальных уравнений первого порядка и систем дифференциальных уравнений первого порядка в Scilab

Для решения дифференциальных уравнений и систем в Scilab предусмотрена функция:

[y,w,iw]=ode([type],y0,x0,x [,rtol [,atol]],f [,jac] [,w,iw])

для которой, обязательными входными параметрами являются: y0 - вектор начальных условий; x0 - начальная точка интервала интегрирования ; x - координаты узлов сетки, в которых происходит поиск решения; f - внешняя функция , определяющая правую часть уравнения или системы уравнений (4.7) или (4.10); y - вектор решений.

Таким образом, для того чтобы решить обыкновенное дифференциальное уравнение вида , необходимо вызвать функцию y=ode(y0,x0,x,f).

Рассмотрим необязательные параметры функции ode:

type - параметр с помощью которого можно выбрать метод решения или тип решаемой задачи, указав одну из строк: adams - применяют при решении дифференциальных уравнений или систем методом прогноза-коррекции Адамса; stiff - указывают при решении жестких задач; rk - используют при решении дифференциальных уравнений или систем методом Рунге-Кутта четвертого порядка; rkf - указывают при выборе пятиэтапного метода Рунге-Кутта четвертого порядка; fix - тот же метод Рунге-Кутта, но с фиксированным шагом; если параметр не указан, то автоматически выбирается либо метод прогноза-коррекции Адамса для нежестких задач, либо вариант “stiff”;

rtol, atol - относительная и абсолютная погрешности вычислений, вектор, размерность которого совпадает с размерностью вектора y, по умолчанию rtol=0.00001, atol=0.0000001, при использовании параметров rkf и fix - rtol=0.001, atol=0.0001;

jac - матрица, представляющая собой якобиан правой части жесткой системы дифференциальных уравнений, задают матрицу в виде внешней функции вида J=jak(x,y);

w, iw - векторы, предназначенные для сохранения информации о параметрах интегрирования, которые применяют для того, чтобы последующие вычисления выполнялись с теми же параметрами.

Рассмотрим использование функции на примере следующих задач.

ЗАДАЧА.

Решить задачу Коши

Решение

Перепишем уравнение следующим образом:

Далее представим его в виде внешней функции и применим функцию z=ode(y0,x0,x,f), в качестве параметров которой будем использовать следующие параметры:

f - предварительно созданная функция f(x, y) (f(x,y)=-y+sin(xy) – правая часть нашего дифуравнения);

x - координаты сетки;

y0, x0 - начальное условие y(0)=1.5;

z - результат работы функции.

Приведем листинг программы, которая строит график, моделирующий процесс, описанный данным уравнением. Здесь х меняется от 0 до 35 с шагом 1. Решение уравнения в точках 0, 1,2,…,35 задается вектором z, координаты которого можно распечатать, убрав запятую в конце 5 строки.

clc

xbasc()

function w=f(x,y),w=-y+sin(x*y),endfunction;

x0=0;y0=1.5;x=0:1:35;

z=ode(y0,x0,x,f);

plot(x,z);xgrid(5)

Графическое решение дифференциального уравнения

Заменив пятую строку на z=ode("rkf",y0,x0,x,f), получим решение уравнения методом Рунге-Кутта.