Текст М-файла

function teplo_1_expm

s=input('s=');

g=input('g=');

n=input('n=');

l=input('l=');

alpha=input('alpha=');

h=l/n;

c=alpha/h^2;

m=n-1; a1=ones(m-1, 1);

A=diag(a1, -1)-2*eye(m)+diag(a1, 1), A=c*A;

u0=g; ul=s;

F=c*[u0; zeros(m-2, 1); ul];

xi=(0: h: l)'; x=xi(2: n);

U0=g+(g+3*s)*x-2*(g+s)*x.^2;

t=[0 0.15 1.5];

nt=length(t); res=zeros(nt, n+1);

fprintf('\n значения функции температуры U(x, t)\n')

for i=1: nt

res(i,:)=[u0 ut(t(i), F, U0, A)' ul];

fprintf('U(%4.2f): ', t(i)), fprintf('%6.2f', res(i,:)), fprintf('\n')

end

hold on

plot(xi, res(1,:), '.-')

plot(xi, res(2,:), 'o-.r')

plot(xi, res(nt,:), '*: g')

grid on

s1=sprintf('t=%2.0f', t(1));

s2=sprintf('t=%4.2f', t(2));

s3=sprintf('t=%4.2f', t(nt));

legend(s1, s2, s3, 0)

title(sprintf('U(x, t)=exp(At)*U0-inv(A)*(E-exp(At))*F\n%s %s %s', s1, s2, s3))

function U=ut(t, F, U0, A)

m=size(A); E=eye(m);

eAt=expm(t*A);

U=eAt*U0-A\(E-eAt)*F;

Замечание. Здесь вычисление функции от матрицы реализуется с использованием стандартной функции expm(At), в которой используется алгоритм

Результаты счета

s=12

g=3

n=8

l=1

alpha=1

A =

-2 1 0 0 0 0 0

1 -2 1 0 0 0 0

0 1 -2 1 0 0 0

0 0 1 -2 1 0 0

0 0 0 1 -2 1 0

0 0 0 0 1 -2 1

0 0 0 0 0 1 -2

значения функции температуры U(x, t)

U(0.00): 3.00 7.41 10.88 13.41 15.00 15.66 15.38 14.16 12.00

U(0.15): 3.00 4.81 6.52 8.03 9.29 10.28 11.02 11.56 12.00

U(1.50): 3.00 4.13 5.25 6.38 7.50 8.63 9.75 10.88 12.00

> >

Дискретно-аналитический метод решения задачи о колебаниях балки при ударе