Казанский Государственный Университет

кафедра вычислительной математики

Выполнил:  студент 3-го курса

Проверил: 

Казань, 2010г

Постановка задачи

Одна из специальных функций математической физики - функция Бесселя, определяется следующим образом

Цель задания - изучить и сравнить различные способы приближенного вычисления этой функции.

Для этого:

1. Протабулировать на отрезке с шагом h с точностью , основываясь на ряде Тейлора, предварительно вычислив его

,

где и получить, таким образом, таблицу

...

...


2. По полученной таблице значений построить интерполяционный полином Ньютона, приближающий


и вычислить погрешность интерполирования

3. На той же сетке узлов построить таблицу приближенных значений используя составную квадратурную формулу трапеций


где а точки разбиения отрезка интегрирования на частей,

Интеграл вычислить с точностью Точность вычисления интеграла определяется сравнением результатов при различном числе разбиения отрезка интегрирования. Именно, точность считается достигнутой, если

4. Построить таблицу обратной к функции

...

...


решая уравнения


следующим итерационным методом: решение уравнения находится по формуле

где 
В качестве начального приближения к точке взять

Выполнение задачи

1) Протабулируем на отрезке с равномерным шагом с точностью , основываясь на ряде Тейлора:

Программная реализация:

Функция J0(x, E) вычисляет значение функции J0(x) по значениям вектора с точностью .

НЕ нашли? Не то? Что вы ищете?

function f = J0(x, E)

f=zeros(size(x));

for i=1:prod(size(x))

  a=1;

  sum=a;

  z=-x(i)^2/4;

  for n=0:1e4

  a=a*z/(n+1)^2;

  sum=sum+a;

  if abs(a)<E, break; end;

  end;

  f(i)=sum;

end;

Использование функции J0(x, E) в main, главном файле программы, для равномерной сетки

a=0;

b=3;

n=10;

eps=1e-6;

x=linspace(a, b,n)

f=J0(x, eps)

% вывод графика функции J0

figure;

plot(x, f,'-r.');

title(['Bessel function ' '; n=' num2str(n)]);

xlabel('x');

ylabel('J_0');

hold on

Результаты:

1.0000

1.3333

1.6667

2.0000

2.3333

2.6667

3.0000

3.3333

3.6667

4.0000

J0

0.7652

0.6026

0.4172

0.2239

0.0376

-0.1276

-0.2601

-0.3514

-0.3972

-0.3971

Рис.1. График функции .

2) По полученной таблице значений построим интерполяционный полином Ньютона, приближающий

Функция NewtonP(x, xi, yi) вычисляет значения полинома Лагранжа на сетке узлов по значениям функции J0(x) в узлах :

%Функция построения интерполяционного полинома Ньютона

function Ln=NewtonP(x, xi, yi)

Ln=zeros(size(x));

Ln(:)=yi(1);

for i=1:prod(size(x))

  P=1;

  for j=1:prod(size(yi))-1

  P=P*(x(i)-xi(j));

  Ln(i)=Ln(i)+P*yi(j+1);

  end;

end;

Использование функции NewtonP(x, xi, yi) в main, главном файле программы, для равномерной сетки:

nn=100;

xx=linspace(a, b,nn);

y=Rr(x, f);

Ln=NewtonP(xx, x,y);

% вычисление погрешности интерполирования err

J0x=J0(xx, eps);

err=norm(J0x-Ln, inf)

error=J0x-Ln;

% вывод графика полинома Ньютона

H=figure;

plot(xx, Ln,'-r');

title(['Newton polinom' '; nn=' num2str(nn)]);

xlabel('x');

ylabel('Ln');

zoom on;

% вывод графика погрешности интерполирования error

ERR=figure;

plot(xx, error,'-r');

title(['ERROR=norm(J_0-Ln)' '; nn=' num2str(nn)]);

xlabel('x');

ylabel('ERR');

zoom on;

Рис.2. График полинома Ньютона.

Рис.3. График погрешности интерполирования.

3) На той же сетке узлов построить таблицу приближенных значений J0(x).

%функция вычисления интеграла

function y=integr(c, d,eps, x)

for i=1:numel(x)

  n=2;

  s=Sn(c, d,n, x(i));

  n=n*2;

  ss=Sn(c, d,n, x(i));

  err=abs(s-ss);

  while err>=eps

  s=ss;

  n=n*2;

  ss=Sn(c, d,n, x(i));

  err=abs(s-ss);

  end;

  y(i)=s;

end;

%функция для вычисления

function s=Sn(c, d,N, xi)

h=(d-c)/N;

s=0;

for i=1:N

  s=s+(podint(xi,(i-1)*h)+podint(xi, i*h))/2;

end; s=s*h;

Использование функции integr(c, d,epss, x)в main, главном файле программы

a=1; b=4;

c=0; d=pi;

n=10;

x=linspace(a, b,n)

epss=1e-4;

J0p=(1/pi)*integr(c, d,epss, x)

Рис.4. График полинома J0(x).

Рис.5. График погрешности интерполирования.

4) Построить таблицу обратной к функции

%Функция для вычисления обратной

function [out, F]=ObrFunc(f0,fn, eps, x)

count=prod(size(x));

z0=x(1);

for i=1:count

  F(i)=f0+i*(fn-f0)/count;

  err=1+eps;

  while (err>=eps)

  k1=J0(z0,eps)-F(i); k2=J0_2(z0,eps);

  z11=-k1/(4*k2^2);

  z12=7*k2-3*J0_2(z0-(2*k1/3*k2),eps);

  z1=z0+z11*z12;

  z0=z1;

  err=abs(J0(z1,eps)-F(i));

  end;

  out(i)=z1;

  z0=z1;

end;

%Функция для нахождения

function f=J0_2(x, E)

f=zeros(size(x));

  sum=-1;

  z=-x^2/4;

  a=-z^(-1);

  for n=2:1e4

  a=a*z*x/2;

  sum=sum+a;

  a=a/(n-1);

  if abs(a)<E, break; end;

  end;

  f=sum;

Использование функции ObrFunc(f0,fn, eps, x)в main, главном файле программы

  a=1; b=4;

  n=10;

  eps=1e-6;

  x=linspace(a, b,n);

  f=J0(x, eps);

  [z, F]=ObrFunc(f(1),f(n),eps, x);

Рис.6. График обратной функции.