Tein sinkert-funktiosta uuden version: siniSkertoimet. (Vanha on ennallaan.)
Tässä on uusin ns05.mpl (Sama on myös /p/edu/mat-1.192/:ssa.)
siniSkertoimet:=proc(f,L,N,tiedosto)
# b:=siniSkertoimet(f,L,N,tiedosto); Funktion f sinisarjan kertoimet.
# b:=siniSkertoimet(f,L,N); palauttaa b-kertoimet, mutta ei kirjoita
# tiedostoon.
# f -- funktio
# L -- jakson puolikas
# N -- Mihin saakka kertoimia lasketaan
# tiedoston nimi merkkijononona (ilman polkua), esim. "bkert.txt"
# Kertoimet palautetaan liukulukutaulukkona (array) ja kirjoitetaan tiedostoon,
# josta ne voi suoraan lukea Matlabiin load-komennolla:
# Lb=load('/tmp/bkert')
# Kirjoitetaan vektori [L,b1,b2,...,bN].
# Työtä voi jatkaa joko Maplessa tai Matlabissa.
Nyt siis kirjoitetaan /tmp/-hakemistoon käyttäjän nimeämään tiedostoon.
>> Lb=load('/tmp/bkert');
>> L=Lb(1)
>> b=Lb(2:end);
Muutin lamposin-funktion sellaiseksi, että se kutsuu täsä uutta.
Siinä siis käyttäjälle näkyvä muutos on, että tiedostonnimi pitää antaa.
>> x=linspace(0,2*pi,30); >> y=x; >> [X,Y]=meshgrid(x,y); >> surf(x,y,Z) % Onhan niitä! >> Z=sin(3*X.*Y); >> M=max(Z(:)) % Jonoutetaan M pitkäksi vektoriksi (sarakkeittain) >> xmaximit=X(Z>M-0.01); ymaximit=Y(Z>M-0.01); % Ideksoidaan ehtomatriisilla >> plot(xmaximit,ymaximit,'o') >> isotarvot=sin(3*xmax.*ymax); % Tarkistus.
> restart:with(plots):setoptions3d(orientation=[-30,50],axes=box):with(LinearAlgebra):
> # read("/p/edu/mat-1.192/ns05.mpl");
> f:=x->...:
# Kehitä sinisarjaksi.
> op(sinkert);
> op(ParitonJ);
> plot(ParitonJ(f),-1..1); # Muistele
> plot(JJ(ParitonJ(f),-1..1),-3..3); # näitä
Muodosta b-kertoimet sinkert-funktiolla (vaikka 50 kerrointa, kestää hiukan).
Laske myös analyyttisesti, kuten ennen. (sinisarjan kertoimet saat suoraan yllä olevasta sinkert-koodista)
> L:=1:omega:=Pi/L:banal:= ... > banal:=trigsiev(banal,n); > ba:=[seq(banal,n=1..10)]; > evalf(%); % Vertaa. > seq(b[j],j=1..10); > osasumma:=unapply(sum(b[n]*sin(n*Pi*t/L),n=1..N),b,t,N): > osasumma(b,x,6); > osasumma(ba,x,6); > osasumma(b,1/2,6); > addversio:=(b,t,N)->add(b[n]*sin(n*Pi*t/L),n=1..N): > addversio(ba,1/2,6); > evalf(%);Huom! Jos haluttaisiin pelata täysin varman päälle, lisättäisiin vielä jakson puolikas L osasumma-funktion argumenttilistaan, mutta nämä ovat yleensä hyvin vuorovaikutteisia "pikamääritysfunktioita". (b:n lisääminen argumenttilistaan ja sen käsitely yo. muodossa tekee hommasta turvallisemman kuin aiempi tapa.)
> display(plot(JJ(ParitonJ(f),-1..1),-3..3),plot(osasumma(b,x,15),x=-3..3,color=blue));
> lprint(banal);
> writeto("/tmp/bf.txt"):
> lprint(banal);
> writeto(terminal):
> banal;
Matlab-ikkuna
clear; close all
x=linspace(-2,2);
[summa,osasummat]=Msinsarja(1,x);
plot(x,summa); figure; plot(x,osasummat)
figure(1); hold on; fplot(inline('x.^2'),[0,1],'r') % Jos funktiomme oli x->x^2
% Kokeillaan 3-argumenttista muotoa, sitä varten luetaan kertoimet tiedostosta:
b=load('/tmp/Fb.txt');
summa=Msinsarja(1,x,b(1:15)); % Otetaan vain 15 ensimmäistä kiinnostutaan vain
% summasta (1. paluuargumentti).
plot(x,summa) % Hauskaa!
f:=x->piecewise(...); plot(f,0..1): # Funktio voidaan piirtää myös näin. Muodostataan ensin analyyttinen ratkaisu (kuten teht. 3:ssa) > S50:=lamposin(...): > f2:=(x,t)->f(x): # Piirretään huvin vuoksi "alkulämpöteltta" ympärille. > display(plot3d(f2(x,t),x=0..1,t=0..0.2,style=wireframe),plot3d(S50(x,t),x=0..1,t=0..0.2));Nyt siihen differenssiin:
> nx:=49: dx:=1/(nx+1): dt:=1/5000: nt:=500: nt*dt;V:= ...; > matrixplot(Matrix(V),style=wireframe); > > kuva:=k->plot([seq([j*dx,V[k][j]],j=1..nx)]); > display(seq(kuva(k),k=1..nt),insequence=true); > x:=[seq(j*dx,j=1..nx)]: x:=evalf(x); > > t1:=evalf(nt*dt); > tarkat:=map(x->S50(x,t1),x); > V[nt]; > tarkat-V[nt]; > max(op(%));Virhe diskretointipisteissä luokkaa ...
> display(plot(zip((a,b)->[a,b],x,tarkat)),plot(zip((a,b)->[a,b],x,V[nt]),color=blue)); Sivumennen sanoen voitaisiin määritellä Matlab-tyylinen plot: > matlabplot:=(xlist,ylist)->plot(zip((a,b)->[a,b],xlist,ylist)): Entä jos r > 1/2 ? Ota vaikka dt:=0.000201: ja tee samat > matlabplot(x,tarkat-V[nt]); > display(matlabplot(x,V[nt]),matlabplot(x,tarkat));
function U=lampodiffB(RAA,L,nx,nt,dt) % Lämpöyhtälön ratkaisu implisiittisellä menetelmällä. % Syötteet: % % RAA -- Reuna-ja alkuarvot, funktio, joka palauttaa kolme tulosta... % Esim: Kirjoita tämä tiedostoon RAA.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % function [vreuna,aa,oreuna]=RAA(x,t) % % Testaa ensin tämä, vaikka: % vreuna=1./t; aa=(sin(pi*x/2)).^2; oreuna=t; % % [v,a,o]=Raa(1:10,-1:1) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % L,nx,nt,dt kuten aina. % % Testaa sitä mukaan kuin kirjoitat rivejä: % L=1; nx=10; dt=0.1; nt=3; % A=lampodiffB(@RAA,L,nx,nt,dt); full(A) % Kohta 1 dx=1/(nx+1); r=dt/dx^2; e=ones(nx,1); A=spdiags([-e 2*e -e],-1:1,nx,nx); U=A % Kohta 1. ....Kaikkea testailua ei kannattane kirjoittaa tähän, vaan esim harj5teht6.m-tiedostoon. Tosin yllä hahmoteltukin voi olla kätevää. Kirjoita istunnossa help lampodiffB ja jatka ...
Suorita yhtälösysteemin ratkaisu LU-hajotelmaa käyttäen, kuten vaikka heat3.m:ssä