Ratkaisut: harj11 LV

Teht. 1

  u(x,t)=5 sin(2x)e-4t- 30sin(3x)e-18t

close all; clear
x=linspace(0,pi,30);
t=linspace(0,0.4,10);
[X,T]=meshgrid(x,t);
U=5*sin(2*X).*exp(-4*T)-30*sin(3*X).*exp(-18*T);
surfc(x,t,U); shg
colorbar
figure
plot(x,U); shg

Teht. 2

> u(x,t)=400/Pi*sum(1/n*sin(n*Pi*x/10)*exp(-lambda[n]^2*t),n=1..infinity); 

                         /infinity     n Pi x                2   \
                         | -----   sin(------) exp(-lambda[n]  t)|
                         |  \            10                      |
                     400 |   )     ------------------------------|
                         |  /                    n               |
                         | -----                                 |
                         \ n = 1                                 /
             u(x, t) = -------------------------------------------
                                           Pi

function [Z,x,T]=osdysin(L,c,b,T,M)
% HA 25.11.03 
% Lasketaan lämpöyhtälön u_t = c^2 u_xx sinisarjaratkaisu korkeusmatriisina
%
% Input: L:   jakso/2
%        c:   c (Lämpöyhtälön c)  
%        b:   sinisarjan kerroinvektori (nollat mukana)
%        T:   aikapisteet
%        M:   x-diskrertointi (oletus 50)
%        
% Output: [Z,x,t] -  lämpötilamatriisi, x-vektori, t-vektori
%
% Esim:
% clear; close all
% L=10;c=sqrt(1.14);nb=10;
% b=1./(1:nb);b=rem((1:nb),2).*b;
% kerroin=400/pi; b=kerroin*b;
% T=linspace(0,10,10);
% [U,x,t]=osdysin(L,c,b,T);
% surfc(x,t,U);  % Pintapiirros
% figure
% plot(x,U)      % Lämpötilaprofiilit t-vektorin 
%                %  ilmaisemilla ajanhetkillä.
%
N=length(T);
if nargin==4, M=50; end
nb=length(b);
x=linspace(0,L,M);
Z=zeros(M,N);
for n=1:nb;
   Z=Z+(b(n)*sin(n*pi*x./L))'*exp(-((n*pi*c/L)^2)*T);
end;
Z=Z';

close all; clear
L=10;c2=1.14;
lam2=c2*pi^2/L^2

% Rakennetaan Fourier-sinisarjan kerroinvektori. 
% Huom! Nollat pitää olla mukana
N=20;
b=1./(1:N);  % Kokonaislukujen käänteisluvut
b(2:2:N)=0;  % Parillisilla indekseillä 0.
b=400/pi*b;
c=sqrt(c2); 
T=0:10;
[U,x,t]=osdysin(L,c,b,T,100);
surf(x,t,U); shg
colorbar
figure
plot(x,U); shg

Tehtävä 3

Tehtava 4

function A=lampomatr(n,r)
% 25.11.03
% Muodostetaan lämpöyhtälön differenssimenetelmämatriisi,
%jolla iteroimalla saadaan sarake kerrallaan
% uusia lämpösarakkeita.
%   Syötteet:
% r = menetelmän parametri, r=c^2*k/h^2.
%     (Stabiilisuusehto: r <= 1/2 .)
% n = osavälien lukumäärä. (pisteitä n+1, sisäpisteitä n-1)
%    Tulos:
% Palautetaan (n-1)x(n-1)-matriisi
%
D=diag((1-2*r)*ones(1,n-1));
Yla=diag(r*ones(1,n-2),1);
Ala=diag(r*ones(1,n-2),-1);
A=Ala+D+Yla;

>> syms r
>> lampomatr(4,r)
 
ans =
 
[ 1-2*r,     r,     0]
[     r, 1-2*r,     r]
[     0,     r, 1-2*r]
 
>> lampomatr(10,r)
 
ans =
 
[ 1-2*r,     r,     0,     0,     0,     0,     0,     0,     0]
[     r, 1-2*r,     r,     0,     0,     0,     0,     0,     0]
[     0,     r, 1-2*r,     r,     0,     0,     0,     0,     0]
[     0,     0,     r, 1-2*r,     r,     0,     0,     0,     0]
[     0,     0,     0,     r, 1-2*r,     r,     0,     0,     0]
[     0,     0,     0,     0,     r, 1-2*r,     r,     0,     0]
[     0,     0,     0,     0,     0,     r, 1-2*r,     r,     0]
[     0,     0,     0,     0,     0,     0,     r, 1-2*r,     r]
[     0,     0,     0,     0,     0,     0,     0,     r, 1-2*r]

Tehtävä 5

function [u,x,t]=lampodiff(u0,tloppu,r,c) 
% Sauva olkoon 1:n pituinen
% Esim: u_t=u_xx, alkuehto: u(x,0)=f(x)=4*x-4*x^2, 0-reunaehdot.
%  c=1;r=0.5;  % n=5, h=0.2
% f=inline('4*x.*(1-x)')
% u0=f(0:0.2:1);
% [U,x,t]=lampodiff(u0,0.2,r,c);
%
n=length(u0)-1;
h=1/n;  % x-askel
k=r*h^2/c^2; % t-askel (pisin, jolla stabiilisuus taattu) 
x=0:h:1; t=0:k:tloppu;
m=length(t);
u=zeros(n-1,m); % alueen sisällä vain n-1 x-jakopistettä
u(:,1)= u0(2:n)'; % jätetään päätyjen nollat pois 
M=lampomatr(n,r);
for j=1:m-1
   u(:,j+1)= M*u(:,j);
end
 % Liitetään lopuksi matriisin ylle ja 
 % alle nollarivi, jotta saadaan reuna-arvotkin mukaan.
nollarivi=zeros(1,m);
u=[nollarivi;u;nollarivi];

Otetaanpa tähän:

>> help lampodiff

  Sauva olkoon 1:n pituinen
  Esim: u_t=u_xx, alkuehto: u(x,0)=f(x)=4*x-4*x^2, 0-reunaehdot.
% Suoeitetaan tästä alkaen:
  c=1;r=0.5;  % n=5, h=0.2
  f=inline('4*x.*(1-x)')
  u0=f(0:0.2:1);
  [U,x,t]=lampodiff(u0,0.2,r,c);
% U:n sarake j kuvaa sauvan lämpätiloja hetkellä t_j.
% Jotta saataisiin oikeassa asennossa oleva taulukko, transponoidaan:
>> U'

ans =

         0    0.6400    0.9600    0.9600    0.6400         0
         0    0.4800    0.8000    0.8000    0.4800         0
         0    0.4000    0.6400    0.6400    0.4000         0
         0    0.3200    0.5200    0.5200    0.3200         0
         0    0.2600    0.4200    0.4200    0.2600         0
         0    0.2100    0.3400    0.3400    0.2100         0
         0    0.1700    0.2750    0.2750    0.1700         0
         0    0.1375    0.2225    0.2225    0.1375         0
         0    0.1113    0.1800    0.1800    0.1112         0
         0    0.0900    0.1456    0.1456    0.0900         0
         0    0.0728    0.1178    0.1178    0.0728         0

% Jos halutaan nähdä eteneminen alhaalta ylöspåin, voidaan tehdä:

flipud(ans)

ans =

         0    0.0728    0.1178    0.1178    0.0728         0
         0    0.0900    0.1456    0.1456    0.0900         0
         0    0.1113    0.1800    0.1800    0.1112         0
         0    0.1375    0.2225    0.2225    0.1375         0
         0    0.1700    0.2750    0.2750    0.1700         0
         0    0.2100    0.3400    0.3400    0.2100         0
         0    0.2600    0.4200    0.4200    0.2600         0
         0    0.3200    0.5200    0.5200    0.3200         0
         0    0.4000    0.6400    0.6400    0.4000         0
         0    0.4800    0.8000    0.8000    0.4800         0
         0    0.6400    0.9600    0.9600    0.6400         0

>> surf(x,t,U'); colorbar
figure
>> plot(x,U')