Forum Mjukvara Programmering och digitalt skapande Tråd

Matlab uppgift: The Bratu Problem

1
Skriv svar
Permalänk
Medlem
●

Matlab uppgift: The Bratu Problem

The Bratu problem is the boundary value problem -u'' = lambda*e^u; 0 < t < 1
for lambda > 0, subject to the boundary conditions, u(0) = u(1) = 0:

Depending on the parameter lambda, this nonlinear problem has two, one or no solutions.
For lambda = 1, it has exactly two solutions. One of them is easily computable,
the other one is hard to find. Your task is to approximate the solutions for lambda = 1
by different numerical methods.

1. Use the finite difference method with n equidistant interior nodes! Set h =
1=(n + 1) and ti = ih, for i = 0; : : : ; n + 1. The resulting system becomes
nonlinear in the discrete approximations yi. Plot the sequence of solutions
which you obtain for n = 1; 3; 7; and 15!

2. Try to find the second solution! Hint: Both solutions have positive initial
slope.

Solution:

m-fil 1:
function f = bratufkn(u);
n =15;
h = 1/(n+1);
A = eye(n)*(-2/h^2);

H = ones(1,n-1)*1/h^2;
A = A+diag(H,-1)+diag(H,1);

f = -A*u-exp(u);

m-fil 2:
options = optimset('TolFun',1e-8,'Display','iter');

n = 15;
A = ones(n,1);
A(n/2+.5,1) = 10

[u,fval,exitflag,output] = fsolve(@bratufkn,A,options)
u = [0;u;0];
h = 1/(n+1);
l = size(u);

for i = 0:n+1
l(i+1) = i*h;
end

plot(l,u,'r-',l,u,'go')
xlabel('t')
ylabel ('u')
title (' n = 15')

Nu är det så att jag har fått m-filerna av en kamrat men skulle behöva hjälp med att förstå de. Någon som vill bolla lite?

Redigera
Citera flera Citera
Permalänk
Medlem
●

Har lite brått men du kan väl börja med att skriva lite vad du själv tänker om uppgifter. Är det en ren matlab-kurs eller är det en kurs i numeriska metoder?
Vet du hur du diskretiserar olika operatorer, så som derivata, andra derivata?
Vet du hur du gör diskretiseringen så att den är "optimal" ur den synpunkten att din numeriska lösning kommer att vara konvergent samt konsistent och vidare stabil?

I din första m-fil så deklareras funktionen bratufkn. Som tar hand om diskretisering av andraderivatan med matrisen A.
En bra övning torde vara att göra denna diskretisering med hjälp av spdiag istället vilket kommer resultera i ett snabbare program då matlab kan nyttja glesa matriser istället för fulla.
följande två rader kanske kan hjälpa dig... (hämtat från dokumentationen till spdiag)
e = ones(n,1);
A = spdiags([e -2*e e], -1:1, n, n);
om du vill jämföra utseende så kika på den med t.ex. full(A)

andra filen får du hjälp med om +24h eller tidigare om någon annan förbarmar sig.

Som sagt, vore trevligt att få se vad du själv gör för arbete, det är ju du som ska kunna detta

Visa signatur

weeeee

Redigera
Citera flera Citera
Permalänk
Medlem
●
Skrivet av mounte:

Har lite brått men du kan väl börja med att skriva lite vad du själv tänker om uppgifter. Är det en ren matlab-kurs eller är det en kurs i numeriska metoder?
Vet du hur du diskretiserar olika operatorer, så som derivata, andra derivata?
Vet du hur du gör diskretiseringen så att den är "optimal" ur den synpunkten att din numeriska lösning kommer att vara konvergent samt konsistent och vidare stabil?

I din första m-fil så deklareras funktionen bratufkn. Som tar hand om diskretisering av andraderivatan med matrisen A.
En bra övning torde vara att göra denna diskretisering med hjälp av spdiag istället vilket kommer resultera i ett snabbare program då matlab kan nyttja glesa matriser istället för fulla.
följande två rader kanske kan hjälpa dig... (hämtat från dokumentationen till spdiag)
e = ones(n,1);
A = spdiags([e -2*e e], -1:1, n, n);
om du vill jämföra utseende så kika på den med t.ex. full(A)

andra filen får du hjälp med om +24h eller tidigare om någon annan förbarmar sig.

Som sagt, vore trevligt att få se vad du själv gör för arbete, det är ju du som ska kunna detta

Gå till inlägget

Roligt att du vill bolla lite.
Jo som du säger så måste jag ju kunna det här till slut så det är bra att jag lär mig och inte bara får svaret. Jag läser en restkurs i numeriska metoder. Jag har egentligen inga vidare kunskaper i Matlab förutom enkla uträkningar, matriser och lite kod som fsolve, function med variabler, och några till. SÅ där ligger jag nu. Det är då givet att jag inte vet hur man diskretiserar olika operatorer och så. Det som står i m-filerna är för det mesta ganska främmande. Men vi kan ta det sakta?

Om vi börjar med första m-filen:

function f = bratufkn(u);
n =15;
h = 1/(n+1);
A = eye(n)*(-2/h^2);

H = ones(1,n-1)*1/h^2;
A = A+diag(H,-1)+diag(H,1);

f = -A*u-exp(u);

Här har vi olika variabler och förutsättningar som är angivna i uppgiften. Men när vi kommer till
A = eye(n)*(-2/h^2) blir det svårt.
Du säger att det är en matris och i funktion browser står det att eye står för identity matrix. Ok. Men varför har man använt just denna matris-baserade ekvation?
H = ones(1,n-1)*1/h^2;
A = A+diag(H,-1)+diag(H,1);
I matlab funktion browser står det att "ones" creates an array of all ones. Jag förstår inte riktigt detta. Sista ekvationen är A igen som var en matris men nu är något annat som jag inte förstår alls.

Jag har egentligen inte mycket kvar av kursen bara en inlämningsuppgift på några problem...som jag låtit hänga allt för länge. I mitt fall tror jag inte att det är värt att lära mig allt kursen går ut på från scratch utan bara det som krävs för att förstå uppgiften. Vad tror du?

Redigera
Citera flera Citera
Permalänk
Medlem
●

Hej igen,
för att inte skriva ihjäl mig så tar jag lite hjälp av wikipedia. Om du kikar på: http://en.wikipedia.org/wiki/Finite_difference
så ser du att man kan diskretisera andraderivatan u'' som:
u(x)'' ~ (u(x+h) - 2u(x) + u(x-h)) /h^2 , där h är "litet"
om vi betraktar u som en diskret vektor t.ex. u=[1 2 3 4 5]
så kan du således approximera andraderivatan för u(2) = 2 enligt (med h=1) u''(2) ~ (u(2+1) - 2u(2) + u(2-1)) / 1^2 = 3 - 2*2 + 1 = 0 ( vilket stämmer bra eftersom det u beskriver en linje).
Om du kikar på koefficienterna framför u(x+h), u(x) ... etc. så för andra derivatan är dom 1, -2, 1, eller som en vektor [1, -2, 1], nu får vi inte missa att allt skulle delas med h^2 således blir vektorn [1, -2, 1]/h^2
Om du nu kikar på din första m-fil så kommer A = eye(n)*(-2/h^2) att skapa koefficienten för -2/h^2.
H = ones(1,n-1)*1/h^2 skapar de två andra koefficienterna och allt slås sedan ihop till matrisen A som nu ser ut något i stil med:

           [-2, 1, 0, 0, ...
           [1, -2, 1, 0,... 
A= 1/h^2 * [0, 1, -2, 1, 0, ...
           [...

Om du skulle multiplicera denna matris med u som är vår funktion, så får du A*u = ... ja lika med vad?
om du tar rad för rad så finner du att första raden i A multiplicerad med u motsvarar andraderivatan för u(1) dvs u''(1),
på samma sätt men med rad 2 från A så blir det u''(2) osv.
Så A är den matris som diskretiserar andraderivatan. Det är en tridiagnoal matris och den är således gles. Matlab kan hantera mycket större glesa matriser än fulla matriser. Därför gav jag tipset att använda spdiags för att skapa denna matris.

Jag återkommer med mer information om ett gäng timmar

Visa signatur

weeeee

Redigera
Citera flera Citera
Permalänk
Medlem
●

Tack. Lärlingen återkommer när han studerat in din text nogrant.

Redigera
Citera flera Citera
Permalänk
Medlem
●

Just det, vill bara slänga in lite ljus i tunneln... matlab fixar mycket av sig självt, t.ex. så kan du lösa ditt problem med den inbygda bvp-solvern med t.ex.

bratuode = @(x,y) [y(2); -exp(y(1))]
bratubc = @(ya, yb) [ya(1); yb(1)]

solinit = bvpinit(linspace(0,1,5),[0.1 0]);
sol1 = bvp4c(bratuode,bratubc,solinit);
solinit = bvpinit(linspace(0,1,5),[3 0]);
sol2 = bvp4c(bratuode,bratubc,solinit);
plot(sol1.x,sol1.y(1,:),sol2.x,sol2.y(1,:))
xlabel('x'); ylabel('y');

för lambda = 1

Visa signatur

weeeee

Redigera
Citera flera Citera
Permalänk
Medlem
●

OK. Du ska ha tack för all hjälp. Jag ska studera in allt det här.

Redigera
Citera flera Citera
1
Skriv svar
Senaste nyheterna
Hårdvara
Mjukvara
Övrigt
Nytt i forumet
Datorkomponenter
System
Ljud, bild och kommunikation
Spel och mjukvara